The use of a priori information in ICA-based techniques for real-time fMRI: an evaluation of static/dynamic and spatial/temporal characteristics

Introduction

Real-time fMRI (rt-fMRI) is an emerging neuroimaging tool based on the estimation of brain activity in real-time (typically around 1–2 s; Weiskopf et al., 2004, 2007; deCharms, 2008; LaConte, 2011). This tool can be used not only for overall monitoring of fMRI data quality (Weiskopf et al., 2007) but also for manipulating the cognitive state of the subject based on their own brain activity (Shibata et al., 2011). The neurofeedback approach has been used in various fields of cognitive neuroscience such as attention (Thompson et al., 2009) and emotion (Posse et al., 2003). Neurofeedback approaches have also been used with rt-fMRI in clinical research, such as the study of control of chronic pain (deCharms et al., 2005) and the control of craving (Chiu et al., 2010; Anderson et al., 2011).

Since its advent, rt-fMRI has had to face a number of technical challenges, mainly due to the computational load of the data analysis which directly competes against the goal of providing real-time feedback (i.e., <1 TR). However, recent technological advancements have provided a way to overcome this issue by making large scale computations possible even on standard platforms (Weiskopf et al., 2007; Weiskopf, 2012). These technical advances have enabled us to shift our focus of attention from technical issues to data analysis aspects.

The usual goal of a real-time system is to permit the identification and monitoring of an activity of interest during its ongoing development and actuation. The identification is defined as an initialization phase where the real-time analysis and derived spatial-temporal features to be monitored are defined, usually with a functional localizer (FL) or a classification training step (LaConte et al., 2007). The monitoring represents the execution of the on-line analysis of the event of interest and the real-time delivery of results that can eventually operate on the stimulation paradigm. From a conceptual point of view, it is thus possible to discriminate the identification and monitoring phases and to develop different algorithms and strategies to deal with them.

The initial and still most common analysis framework for rt-fMRI is based on univariate hypothesis-driven approaches, with adaptation of standard algorithms, such as the general linear model family (GLM), to the real-time domain (Cox et al., 1995; Gembris et al., 2000; Hinds et al., 2011). These methods are common mostly because they are associated with ease of interpretability and fast computation. In these approaches both the identification and monitoring phases are typically implemented using hemodynamic response-based models of the expected cognitive tasks and eventual nuisance variables taking place during the rt-fMRI experiment.

Another family of data analysis techniques is represented by the multivariate data-driven algorithms, which have shown a great capability of exploiting the full information content intrinsically present in the data to be analyzed without assuming the explicit shape or timing of the hemodynamic response to a stimulus (McKeown et al., 1998; Mouro-Miranda et al., 2005; Norman et al., 2006). The driving motivation behind these methods is that they allow characterizing functions that may not be detectable without exploiting both second order (variance) and higher-order statistics, thus relying on a greater amount of information. These properties make the multivariate data-driven techniques very appealing for use in the real-time domain. Within this concept several machine learning algorithms have been successfully adopted and exploited in the real-time data analysis framework. The most successful implementations are based on support vector machines (SVM; LaConte et al., 2007; Magland et al., 2011; Sitaram et al., 2011). SVM provides a powerful solution to a number of applications that are subject specific, at the cost of training the classifier and imposing some interpretability issues on the results. In this context, the two phases of the canonical rt-fMRI framework are represented by the two steps of a classifier, i.e., the first phase is the training of the classifier, and the second phase is the test or execution of the classifier (i.e., the classification itself).

In addition to SVM, independent component analysis (ICA), another multivariate data-driven technique, has proven to be very effective in fully exploiting the complete amount of information which is present in the data. ICA enables the extraction of knowledge other than that merely modeled in a classical univariate approach (Hyvrinen and Oja, 2000; Calhoun et al., 2001; Beckmann and Smith, 2004). Furthermore, ICA methods can also be applied in a series of problems for which univariate inference cannot offer a solution, i.e., in experiments that lack a regressor model to be adopted in the univariate analysis. This is the case for resting data analysis or also experiments with particular patient populations (Calhoun et al., 2009).

The idea of translating ICA properties to a real-time implementation was firstly proposed by Esposito et al. (2003) in a seminal paper and implemented as a plug-in in Turbo Brain Voyager software (Goebel, 2012). In this initial work the authors presented a FastICA based rt-fMRI analysis tool exploiting precise design choices and including an identification phase and a monitoring phase. The first identification phase solved the problem of ranking ICs of interest, i.e., a canonical univariate functional localizer step was implemented to define areas of interest. Other ways to solve the problem of ICs ranking could be represented by exploitation of expected characteristic features of the ICs of interest via a classifier (DeMartino et al., 2007). The second monitoring phase used on-line execution of FastICA (implemented in a sliding window fashion) for extracting different ICs. The ICs were ordered on the basis of their spatial overlap with the IC of interest, which in this case consisted of single-slice representation of motor activity derived from a finger tapping localizer.

The work presented by Esposito et al. (2003) was recently extended to evaluate the performance of 14 different ICA algorithms considering as additional variables the model order and different types of a priori knowledge (spatial/temporal; Soldati et al., 2013). This work showed that ICA algorithms such as EVD, amuse, jadeopac, and FastICA were suitable when implemented in the identification phase via a functional localizer since they performed well even without extensive use of a priori knowledge. It is interesting to note that FastICA algorithm represented a good trade-off and its performance was valid in both functional localizer and dynamic monitoring phases. Other algorithms like constrained ICA performed worse without a priori knowledge and may thus be more suited for the dynamic monitoring phase due to their ability to incorporate a priori knowledge. Such a priori knowledge may help guiding the algorithm to detect a specific target IC with higher priority over the other ICs present in the data. However, there are several types of prior information that are available including spatial domain, the temporal domain, or both, and any of these could be used in different ways (as constant references from a localizer or derived dynamically). It is however not clear how these various ways of using priors may affect the performance of the results both in terms of computation time and correlation to a reference optimal ICA. Moreover, the ICA algorithm (FastICA) is stochastic, which means that multiple repetitions of the analysis on the same dataset can give slightly different results, both in the spatial and temporal domains. The problem has been extensively discussed in the literature, with one of the main proposed solutions being based on multiple ICA runs and clustering of the obtained components, with the aim of reducing the issue of stochastic variability (Himberg et al., 2004). Such instabilities can be characterized by the standard deviation of the derived (STD) results (spatial and/or temporal) when the analysis is repeated multiple times on the same dataset. The STD can be considered as a stability performance parameter of the algorithm, lower STD algorithms corresponding to more stable ones. This parameter may be particularly relevant if different ICA-based algorithms are to be considered and compared for real-time fMRI, where the analysis is repeated dynamically during data acquisition.

This study extends previous work (Esposito et al., 2003) in two ways. Firstly, the target IC to be monitored dynamically is identified from a functional localizer using an ICA-based method instead of using a GLM of the hemodynamic response. This approach allows the full analysis pipeline to be multivariate and data-driven. Secondly, novel ICA-based algorithms are proposed that introduce different types of a priori knowledge for the dynamic monitoring of ongoing fMRI activity. The main goal of this study was to evaluate how these algorithms perform with respect to an off-line ICA analysis after the acquisition is complete. The a priori information considered was either temporal, spatial, or both spatial and temporal. In addition, the a priori information was considered both in its static version when derived from the functional localizer, as well as dynamic when estimated recursively as the sliding window progresses over the time course throughout the run. The different ICA-based analysis methods proposed here were tested by artificially simulating a real-time fMRI experiment using real fMRI data from a visual motor study (Calhoun et al., 2001). The original data used is unrelated to a real-time fMRI experiment, but was adopted because it is public and offers robust functional activation in well-known anatomical areas. The measures of performance to compare the various real-time methods were based on the following three metrics: (i) spatial and/or temporal correlation between the independent component (IC) estimated dynamically and the target IC derived from the localizer, (ii) computation time, and (iii) intrinsic stochastic variability of the algorithms as estimated from multiple analysis runs.

Materials and Methods

Dataset

One of the aims of this work was to test a variety of ICA implementations in a fashion which can be directly applied to real world conditions. For the simulation of the rt-fMRI experiment we used a dataset coming from a real publicly available fMRI experiment, with tasks that show robust activation in well-known brain networks. We chose to use the data that comes as part of the GIFT package (Calhoun and Adali, 2006) because the ICA characterization of the task-induced activation networks was extensively tested. The dataset is thus publicly available and in the release it is stated that The Johns Hopkins Institutional Review Board approved the protocol and all participants provided written informed consent.

Imaging parameters

Scans were acquired on a Philips NT 1.5-T scanner. A sagittal localizer scan was performed first, followed by a T1-weighted anatomic scan [repeat time (TR) = 500 ms, echo time (TE) = 30 ms, field of view = 24 cm, matrix = 256 × 256, slice thickness = 5 mm, gap = 0.5 mm] consisting of 18 slices through the entire brain including most of the cerebellum. Next, we acquired functional scans over the same 18 slices consisting of a single-shot, echoplanar scan (TR = 1 s, TE = 39 ms, field of view = 24 cm, matrix = 64 × 64, slice thickness = 5 mm, gap = 0.5 mm, flip angle = 90°) obtained consistently over a 3-min, 40-s period for a total of 220 scans. Ten “dummy” scans were performed at the beginning to allow for longitudinal equilibrium, after which the paradigm was automatically triggered to start by the scanner.

Experiment setup

The GIFT package contains three subjects example data-sets that employ a visuo-motor paradigm derived from other studies (Calhoun et al., 2001). The paradigm contains two identical but spatially offset, periodic, visual stimuli, shifted by 20 s from one another. The stimuli consisted of an 8 Hz reversing checker-board pattern presented for 15 s in the right visual hemi-field, followed by 5 s of a central asterisk fixation, followed by 15 s of checker-board presented to the left visual hemi-field, followed by 20 s of a central asterisk fixation. The 55 s set of events was repeated four times for a total of 220 s. The motor stimuli consisted of participants touching their thumb to each of their four fingers sequentially, back and forth, at a self-paced rate using the hand on the same side on which the visual stimulus is presented.

Pre-processing

The images were first corrected for timing differences between the slices using windowed Fourier interpolation to minimize the dependence upon the reference slice chosen. Next, the data were imported into the statistical parametric mapping software package, SPM99. Data were motion corrected, spatially smoothed with a 6 mm × 6 mm × 10 mm Gaussian kernel, and spatially normalized into the standard Montreal Neurologic Institute space. The data were slightly subsampled to 3 mm × 3 mm × 5 mm, resulting in 53 × 63 × 28 voxels.

In this study the pre-processing steps were not included as part of the real-time fMRI simulations for several reasons: (i) the pre-processed and not the raw data are publicly available as part of the GIFT package (Calhoun et al., 2001) thereby being a reference starting point for various analysis tools, (ii) these pre-processing steps can be performed in real-time as several review studies describe (LaConte, 2011; Caria et al., 2012; Maclaren et al., 2013), (iii) the focus of this simulation study was on the data-driven network characterization through various real-time algorithms. For these reasons, and to keep a manageable number of variables in this study we limit our simulations to the manipulation of real-time analyses that follow the standard pre-processing steps.

Toolbox and PC

The entire simulation work was based on an in-house implementation with MATLAB (2010) of the tested algorithm based on the code of GIFT toolbox (Calhoun and Adali, 2006). Given the ICA algorithms code present in the toolbox, all the data analysis steps (presented in Materials and Methods section) were implemented in an automatic fashion to permit a testing routine to be run by varying parameters, techniques, a priori knowledge, and different subjects. The PC adopted to run the simulations was an Intel(R) Core(TM) i5 CPU M460 @ 2.53 GHz equipped with 6 GB of RAM and running a Windows 7 64-bit OS.

ICA Mathematical Preliminaries

Since all the methods share a common core based on ICA principles, we briefly recall the main concept associated with the ICA. Let’s assume that we have a set of fMRI measurements Y_i where i = 1, …, v is the index of voxels and each Y_i is a vector of y_ij elements, where j = 1, …, t is the index of time points. The entire dataset can thus be represented as a matrix Y of dimensions time points by voxels. Now, let’s assume that the signal measured in the dataset is generated by a subset of n underlying sources which are linearly mixed and summed up. This reflects in the following canonical formulation using the vector-matrix notation.

$Y=AX$ (1)

where Y is the acquired data matrix of dimension equal to the number of time points by the number of voxels, A is the mixing matrix of dimension equal to the number of time points by the number of sources to be recovered and X is the matrix of the sources (i.e., ICs) of dimension the number of sources by the number of voxels. Each jth row of Y is a vector y_ji=1:v representing an fMRI volume in a jth time point and is thus obtained by the linear weighted combination of hidden sources spatial maps y_j = a_j1 x₁ +…+ a_jnx_n∀_j. This means

$Y=\sum _{j=1}^n{a}_j{x}_j$ (2)

Given this definition and assuming that the sources x_n are mutually independent, it is possible to recover those hidden sources by computing an estimation of the unmixing matrix W = A⁻¹ such that

$X=WY$ (3)

is an estimate of the sources. The estimation of W can be obtained via different algorithms, leading to different ICA implementations with different properties and effectiveness (see Bell and Sejnowski, 1995 for details). In this paper the selection of FastICA (Hyvrinen and Oja, 2000) as the core ICA algorithm has been driven by a recent study that compared the performance of 14 different ICA algorithms, and found FastICA to be amongst the most stable ones (Soldati et al., 2013). The FastICA algorithm exploits the non-Gaussianity as a metric of independence of the sources. This means, in the simplified iterative algorithm for several units, that the estimation of W is obtained through the following steps

1. initialize randomly W

2. given $W=\frac{W}{\sqrt{∥W{W}^T∥}}$

3. repeat until convergence $W=\frac{2}{3}W-\frac{1}{2}W{W}^{T}W$ and step 1–3.

Other approximations of the solution can be obtained, but a detailed description of the methods to obtain FastICA decomposition is anyway beyond the scope of the present paper.

Analysis Framework

The purpose of our analysis was to perform an extensive comparison between the standard off-line ICA analysis and several novel on-line ICA methods. The main goals of this study were to evaluate the feasibility of on-line ICA and identify the best performing algorithm from those proposed. For the purpose of our rt-fMRI simulations we proceeded with three different stages: calculation of reference ICs for performance evaluation, calculation of the target ICs from a functional localizer, and estimation of a dynamic IC on real-time. The first stage used FastICA on the complete fMRI time series to identify spatial and temporal IC templates from the networks that were later monitored dynamically. These templates were derived from the full dataset so they were in this sense considered as gold standard references against which the dynamically extracted components were later compared for spatial-temporal accuracy evaluations. The second and third stages were more strictly related to the rt-fMRI simulations. The second stage simulated a functional localizer (FL) session by taking the first 60 TRs of the fMRI time series. FastICA was used on the simulated FL to extract target ICs that was later monitored dynamically. The third stage represented the real-time fMRI simulation, the on-line ICA decomposition that used the information coming from the simulated FL. This last stage of real-time ICA decomposition was performed using the different novel techniques proposed and described in the next subsections.

The proposed framework for performing rt-fMRI used a multivariate and data-driven approach schematically presented in Figure 1. The general structure and workflow can be outlined as follows: (1) The MR data acquired by the scanner was stored during acquisition and made available to the data analysis system as soon as the images were reconstructed. (2) At the beginning of the experiment a short period (typically about 5 min or less) was devoted to acquire data from a FL. In the proposed framework the FL data was analyzed using a blind (unconstrained) ICA algorithm to preserve the multivariate data-driven advantages. Others used univariate methods at this stage (Esposito et al., 2003). (3) An IC of interest was selected from the FL analysis, this IC became the data-informed multivariate ROI whose activity was meant to be monitored dynamically. (4) The IC of interest, along with possible a priori information, could be incorporated in the rt-ICA data analysis algorithm. The ICA algorithm used a sliding window approach and a blind source extraction (BSE) perspective to deliver results at each TR while updating the best match to the target component. This monitored component or other a priori knowledge was then provided recursively to the algorithm, which extracted the actual version of the monitored IC updated by the actual values of data. (5) The monitored IC could be used as in classical rt-fMRI paradigms for visualization, neurofeedback, or brain computer interfaces. The component selected in real-time was the one that has a spatial map which maximally correlates with the reference spatial map component identified during the Functional Localizer (FL) step. The spatial FL component corresponds in turn to the FL component whose temporal correlation with the timing of the paradigm was highest.

FIGURE 1

Figure 1. The experimental design framework. In this structure it is possible to identify two main phases (or steps). The first step is to identify what to monitor, i.e., performing a Functional Localization. This can be done with or without incorporating some kind of a priori knowledge. The second step is to monitor the phenomenon we identified in the previous step using a suitable on-line analysis method.

In this study two main different approaches were investigated to dynamically extract in real-time a target IC: static methods based on Back-projection and dynamic methods based on iteratively performed ICA. In particular, for the dynamic methods, FastICA, and constrained ICA were updated for a sliding window real-time fMRI implementation. The size of the sliding window was fixed and it was chosen based on previous work that systematically evaluated the performance of multiple ICA algorithms as function of window length amongst other variables (Soldati et al., 2013). This study showed that the sliding window length that gave the optimal trade-off between computational speed and spatial/temporal correlation with the results from the whole time course was approximately equal to the period of the behavioral task to be monitored. In our experiment this could be approximated to around 30 s, i.e., 15 TRs. It is worth noting that the size of this window, as pointed out in the discussion, was strongly related to the period of the behavior to be monitored.

Template creation and accuracy estimation

To estimate the accuracy of one technique in correctly describing a monitored IC at one arbitrary time point we generated task related network templates which represented the principal spatial and temporal characteristics of the ICs to be monitored during the simulation. These templates of task related ICs were thus taken into account as reference data to evaluate the quality of rt techniques. This evaluation was obtained via comparing the dynamically reconstructed ICs with these templates using temporal correlation and spatial overlap. To create the templates an ICA analysis was performed on the single subject level by considering all time points (i.e., 220 TR), but using FastICA with the same model order to be used in the on-line implementation (i.e., 5). Three different target ICs were manually selected to simulate their dynamic monitoring see Figure 2: two task related components (RVMT and LVMT) and the task-induced default mode network (DMN).

FIGURE 2

Figure 2. The monitored ICs. An illustrative example of the monitored ICs is reported. Spatial maps of ICs considered in the simulation are obtained from Group ICA 20 ICs. For ease of visualization only the relevant slices are reported here. First column depicts Default Mode Network (DMN). Second and third columns depict the two task related ICs, Right Visuo-Motor Task (RVMT) and Left Visuo-Motor Task (LVMT).

Functional localizer

In a rt-fMRI experiment the functional localizer could be used to identify the IC to be later dynamically monitored. For the analysis of the FL we considered the use of whole brain ICA to maintain a multivariate data-driven method. Previous simulations suggested that the ICA analysis of the FL data is most accurate when using algorithms such as evd, jadeopac, or FastICA (Soldati et al., 2013). We thus performed the FL applying the FastICA algorithm with a model order of 5, the same algorithm used in the template creation, to the first 60 TRs of the time series. With the application of FastICA as FL an unmixing matrix W was estimated of dimension 5 by number of time points. Each row of the matrix represented a time course of a hidden source, and the associated row of the X derived matrix represented the corresponding spatial map. The target IC to be monitored was then automatically selected as the one whose spatial map maximally correlated with the reference template and was monitored later dynamically with the on-line techniques. In this study three components were extracted from the FL for separate evaluations in the dynamic monitoring: the default mode network, the right, and left hemisphere visual motor networks activated by the cognitive task.

On-line techniques

In this section we present the main developed work, that is the methods implemented to perform the on-line monitoring of the sources. To simulate the on-line ICA analysis, coming after the FL, the rest of the time course (i.e., 220–60 TRs) was used to dynamically monitor the FL-derived target ICs using a sliding window approach. Given that the target was to properly exploit the a priori knowledge, different approaches to combine this knowledge and the ICA algorithms were developed. In a comprehensive perspective all the possible combinations were explored. Starting from the concept of sliding window ICA as it was presented by Esposito et al. (2003), more sophisticated and different methods were implemented. The target was to obtain an actual temporal value of the activation of interest and/or an actual spatial map of this component. Two main criteria were the guidelines in these implementations, that is the dynamic of the data and the type of a priori knowledge. The dynamic criterion means how much novel information is exploited and weighted into the on-line method, while the type of information exploited denote the nature of the a priori knowledge, i.e., temporal, spatial, or both. The implementation exploited state of the art ICA algorithms (FastICA, Constrained ICA) with the target of making the implementation easy to reply and distribute.

The following subsections present the details of the different on-line monitoring techniques proposed.

Static method: back-projection. The basic assumption behind this static method was that the brain activation of interest maintains its basic characteristics, in particular its spatial map (SM), relatively stable during the fMRI experiment. If this assumption holds, the spatial ICA performed in the FL step is enough to extract a precise representation of that spatial map that will be later tracked dynamically during the real-time experiment. ICA would be in fact able to create a space described by the directions of the extracted ICs that is fairly representative of the brain state during the performance of the task of interest. Given this assumption, the SM of an IC of interest obtained from the FL can be kept fixed and it should then be possible to simply back project each newly acquired volume of data into this space (i.e., onto the SM of the IC of interest) to be able to quantify the contribution of the new data to the brain activity of interest. This means that no ICA analysis must be performed in real-time, and the results will only depend on the ICA performed in the Functional Localizer session. This contribution will thus represent the time course of the IC. In more detail, the processing steps can be outlined as follows:

1. FastICA was used on the FL to estimate an unmixing matrix W = A⁻¹ and thus the associated SM of the sources X_nIC,v with nIC equal to the number of extracted components (5 in this case) and v equal to the number of voxels.

2. The SM of the desired component was then chosen as the source whose associated time course was the most correlated to the task of the FL, that is we have X_nIC=sel,v . The chosen spatial map was therefore an independent component computed by the FastICA algorithm, which gives it unitary variance and null mean value characteristics.

3. At this point, for each newly acquired volume Y_ith of dimension one by number of voxels we could compute:

$a_{ith}=Y_{i^{th}}{X}_{nIC=sel,v}^{†}$ (4)

where a_ith is the actual single time value of the IC of interest and † denotes the pseudo-inverse. It is worth noting that this can be straightforwardly extended to cases in which multiple components are monitored simultaneously via parallelizing equation (4) for different SM or obtaining a meta-SM via combining different SM, i.e., X matrices.

Dynamic method: recursive temporally constrained. This algorithm is a direct extension of that used by Esposito et al. (2003). The main differences are that here it was applied to the whole brain and that the computation of a priori temporal knowledge was not model-driven, but it was rather data-driven and obtained with an approach based on the previously presented Back-projection method. The actual difference with the previously presented Back-projection method, in which the SM was static, was that we obtained an actual updated dynamic SM via iterative ICA computation. The details of the method are as follows:

1. From the FL a SM was obtained, which was used as in Back-projection method to obtain temporal a priori information in subsequent steps

2. During the experiment

(a) using the Back-projection the time course of the brain activation was extracted

(b) FastICA algorithm with model order 5 and time window length 15 TR was applied to the data with a sliding window approach. The FastICA was temporally constrained using the a priori temporal constraint (obtained using the Back-projection) to initialize the mixing matrix A.

In practice a sliding window of dimension Δ was updated for each newly acquired volume n leading to a matrix Y_[n-Δ,n],v of dimension Δ by number of voxels. This matrix and the SM obtained in the FL step (i.e., X_nIC=sel,v) were used to extract a time course in a data-driven way in the same fashion as for the Back-projection algorithm, resulting in a time course a_n of dimension nTP by one. With the actual data matrix Y_{[n-nTP,n+nTP],v} and the time course a_n it was then possible to apply FastICA to extract the actual SM of the component of interest. This was done by initializing the first entry of the W matrix with the inverse of a_n, given that W = A⁻¹, in the routine presented in the ICA mathematical preliminaries section. The result was the actual SM IC (i.e., Xnew) present in the data whose behavior was closest to the reference time course. In other words the extracted IC was constrained to be as close as possible to the reference one at the initial step, permitting a much more dynamic computation of the IC and thus update of the monitoring. In this approach the SM was dynamically updated each time a new volume was acquired.

Dynamic method: recursive spatially constrained. As in the previous method, also in this dynamic method, the on-line monitoring required a continuous update of the ICA decomposition matrix. There were two main differences with respect to the RTC method: (i) the ICA algorithm was a spatially constrained ICA (Lin et al., 2010), and (ii) the a priori knowledge was spatial instead of temporal. In this approach, the knowledge of an a priori SM of the IC of interest (obtained by the FL) permitted constraining the computation of the ICA algorithm. The constrained ICA algorithm was applied on time windows of data still of length 15 TR with a sliding window approach. The extracted IC, although based on newly acquired data, was forced to be spatially as close as possible to the spatial a priori given map (i.e., to the SM obtained during FL). This means that the dynamically extracted IC represented the SM of the brain activity of interest in the shape that was actually present in the novel data, thus dynamically updated. The associated time value was given by an approach similar to the Back-projection method but depending on the dynamically updated SM, i.e., given the new SM (i.e., Xnew), by computing $a_{ith}=Y_{ith}{X}_{new}^{†}$.

Dynamic method: recursive spatio-temporal method. This algorithm implemented the possibility of obtaining actual dynamic values from both the time course and the spatial map with two concatenated steps. This was obtained by combining the previous methods to obtain a fully updated on-line method based on the following steps:

1. Back project the actual data on the SM of FL, obtaining the actual value of time course in the FL space $a_{FL}=Y_{ith}{X}_{nIC=sel,v}^{†}$ (note that in the BP algorithm we assumed little or no difference between the template space and the actual subject space)

2. Apply the temporally constrained algorithm (i.e., initialize the W matrix exploiting the a_FL) to obtain the SM (i.e., X_sub) in the actual subject space (note that this is different from the FL space)

3. Apply the spatially constrained ICA with the SM in the actual subject space to obtain the actual time course value in the subject space, that is $a_{sub}=Y_{ith}{X}_{sub}^{†}$.

Adopting the described steps it was thus possible to obtain temporal and spatial values of the brain activation of interest fully exploiting the actual data, thus adapting to the dynamic changes which could occur, but keeping as a target the characteristics defined in the FL session.

Variability Effects from the Stochastic Nature of ICA

ICA methods (with some exceptions like the jade algorithm) are typically non-deterministic since there is a stochastic component in the analysis. This introduces variability each time the algorithm is run, which in turn can affect the computation time and the performance of the dynamic monitoring of a target IC. Such variability effects were investigated by repeating the analysis 10 times for each subject on the same data, and then computing the standard deviation across repeated trials for the mean correlation between dynamic and template spatial maps and temporal time course.

Results

Using publicly available fMRI data from a previous experiment (Calhoun et al., 2001) we simulated a real-time acquisition in a sliding window approach to evaluate the performance of four implementations of ICA with different uses of a priori information: (i) Back-projection of constant spatial information derived from a functional localizer (BP), (ii) dynamic use of temporal (RTC), (iii) spatial (RSC), or (iv) spatio-temporal ICA constrained data (RSTC).

Given the stochastic nature of the ICA algorithms used, the variability of the spatial and temporal results was evaluated for each subject on each of the target networks (Default Mode Network (DMN), Right Visual Motor and Left Visual Motor Task related components (RVMT and LVMT respectively)) and for each of the four ICA implementations. The results showed in Figure 3 point out that stochastic effects can introduce variability in the performance of the IC order ranking accuracy up to 10%, sometimes producing large fluctuations. This behavior suggested that none of the four ICA implementations gave consistently the lowest sensitivity to fluctuations due to stochastic effects, although the BP method tended to be the lowest in 15 out of 18 cases.

FIGURE 3

Figure 3. Variability of dynamic tracking performance results due to the stochastic nature of ICA. Performance is here represented by two metrics: spatial or temporal correlation between the template and the dynamically tracked IC, averaged along the time course. The calculations were repeated 10 times for each subject, for each of the three networks evaluated (default mode or DMN, right visual motor or RVMT, and left visual motor or LVMT), and for each ICA implementation (Back-projection or BP, temporally constrained or TC, spatially constrained or SC, spatio-temporal constrained or STC). The variability of the dynamic tracking performance results is expressed as the standard deviation across trials, per subject, brain network, and ICA implementation.

Proceeding with further analysis it was possible to focus on the evaluation of stability of results across subjects and across different monitored ICs, as presented in Figure 4. This figure reports the standard deviation across subjects of the mean (across trials, for each subject) spatial and temporal correlation for the monitored ICs.

FIGURE 4

Figure 4. Variability of dynamic tracking performance results due to subjects. Similar to Figure 2, but here mean results across trials are used to compute variability across subjects, expressed as standard deviation. The subject variability is shown for each of the three networks evaluated (default mode or DMN, right visual motor or RVMT, and left visual motor or LVMT) and for each of the four ICA implementations (Back-projection or BP, temporally constrained or TC, spatially constrained or SC, spatio-temporal constrained or STC).

Finally the performance evaluation numbers were reported in Figure 5 and Table 1, in which one can see the spatial and temporal correlation between the reconstructed time courses and spatial maps of the monitored ICs and the reference templates of those ICs.

FIGURE 5

Figure 5. Overall dynamic tracking performance in reconstructing ICs. The mean and standard deviation performance results are shown across subjects and trials in terms of both spatial and temporal correlation with the template ICs. The results are shown for each of the three networks evaluated (default mode or DMN, right visual motor or RVMT, and left visual motor or LVMT) and for each of the four ICA implementations (Back-projection or BP, temporally constrained or RTC, spatially constrained or RSC, and spatio-temporal constrained or RSTC).

TABLE 1

Table 1. Performance.

The table reports the performances in terms of computational time necessary to update the actual value of time course or spatial maps for each new available data volume. This represents another critical issue of a real-time analysis.

The results show that the Back-projection method offered the highest performance both in terms of time course reconstruction (correlation value to the template time course was significant and quite high, around 0.9), and speed (computation of update value was far below the TR). This method was very fast and effective as long as the monitored IC had a strong and well defined behavior and/or it was well extracted in the FL, since it relied on an accurate description of the spatial behavior. The fluctuation reported in the figures represents error fluctuations in the FL phase which directly reflect in the Back-projection method. The dynamic methods offered comparable performances at cost of higher computational time (CT) (around 2 s for RTC). In particular the spatio-temporal method performed comparably in terms of CT to Back-projection, offering more variable performances in terms of reconstruction of spatial maps and time courses.

Discussion

In the present work we presented and evaluated different methods to combine ICA-based algorithms for real-time fMRI. The motivation for this work was to investigate how the advantages of such multivariate data-driven based methods can be adapted to real-time fMRI applications, extending previous work (Esposito et al., 2003). One goal of this work was to simulate a realistic scenario fully based on ICA consisting of two essential steps. The first step was dedicated to identifying brain networks of interest from the ICA of a functional localizer. The second step consisted in dynamically monitoring a target IC (derived from the first step) with the use of different types of a priori knowledge in the computations. The a priori information considered ranged from static to dynamic, where spatial maps and time courses can be updated separately or together to give more weight to the dynamic monitoring of data within a pre-established time window in the fMRI time course. The incorporation of a priori information was motivated to address the challenge of identifying and keeping track of a specific IC of interest, despite all the other ICs that might be present in the data. This work therefore focused on evaluating different ways of using prior information about the target IC to monitor such that during the dynamic monitoring phase the target IC could be effectively detected with higher priority relatively to other possible ICs.

The ICA-based techniques presented for the on-line monitoring were characterized by different advantages and disadvantages. Overall findings confirmed two general expected features: (i) the dynamic monitoring performance was directly related to the strength of activation of the target IC identified in the functional localizer, stressing the importance of this first step, and (ii) as algorithms became more adaptive in the use of spatial and/or temporal priors in the dynamic monitoring, they introduced less stability in the performance results compared to off-line results. This reflected the intrinsic differences between static off-line analysis and dynamic one.

Back-projection is the only method presented for which the ICA is computed only once in the FL session, and not updated later during the on-line monitoring. This means that this method as implemented here is based on ICA since it depends on the quality of ICA performed in the FL session, but it is not a fundamentally ICA method. Back-projection could in principle also be used by defining a target brain network from the functional localizer with a standard general linear model that makes assumptions on the hemodynamic responses. The use of ICA, however, allowed the analysis to be fully data-driven (Esposito et al., 2003; Beckmann and Smith, 2004; Norman et al., 2006; Calhoun et al., 2009; Magland et al., 2011; Sitaram et al., 2011) and this represents an advantage in all those experimental designs where the classical ICA showed to be robust, as in all cases lacking a defined regressor. The Back-projection technique had the positive features of being stable in terms of lowest fluctuation across trial and subjects, very fast relative to the TR of fMRI data acquisition (since it just involves a matrix multiplication), conceptually simple, and being able to monitor more than one IC of interest. The main potential disadvantage of the Back-projection method was related to its non-adaptivity, since it assumed that the target IC of interest was always present with the same properties, i.e., a fixed spatial map was considered.

The temporally constrained ICA was more adaptive to data with respect to Back-projection. Even if similar to what was presented by Esposito et al. (2003) this method offered different characteristics. The main one was that the reference time course used as constraint was not obtained using a hemodynamic model, but it was extracted from the data in a multivariate data-driven way by the FL. Moreover in this method the reference time course was updated in a similar way to Back-projection, while the crucial difference was that the spatial map updates iteratively each time new data become available. A characteristic of the temporally constrained ICA was that the dynamic spatial map generated was derived from the time course used to initialize the ICA algorithm. This time course, being derived from the Back-projection of actual data on a static space (i.e., keeping the spatial map of the IC of interest fixed), was strictly related to the quality of the FL. For this reason the time course reconstructed was in the template space (i.e., FL space), while the spatial map was in the subject space, being obtained by exploiting the reconstructed time course as a priori knowledge during the application of the ICA algorithm. A limitation of the temporally constrained algorithm was that its mean computation time was more than one order of magnitude higher relative to all the other methods tested. This is due to the fact that the FastICA algorithm adopted in it, while performing generally lower than a TR, sometimes (around 2–3% of the times) got stuck in a local minimum thus increasing the time to perform the decomposition (in some cases from 1.5 upto 8 s). A possible solution to this would be to skip the updating of the information for those volumes which exceed a pre-determined temporal limit to update.

The spatially constrained ICA assumed a fixed spatial map of the IC of interest. This approach suffered from the small amount of data available for the decomposition. The main advantages included low computational time and low variability of the results, qualities that make it a good candidate for use in real-time experiments.

The combined implementation of spatial and temporal constrained ICA permitted a better description of the actual dynamic behavior of data, thus focusing on data characteristics which were strongly transient and for this reason probably not modeled in the off-line static analysis, which privileged extraction of static periodic or quasi-periodic behaviors. This method enabled us to obtain valuable results both in terms of accuracy and computational time. Its main disadvantage was that it was less able to characterize static aspects of the data.

A further consideration is needed related to the variability of monitoring performance. Three kinds of variability were investigated in the simulations. The first one was due to the stochastic nature of principal ICA algorithms, which caused different results to be obtained in different runs of the algorithm on the same data. Multiple repetitions of the analysis showed that this variability can affect computation time, but the obtained performance had a stability better than 10%. The second kind of variability identified was subject specific which caused about 20% of the variability. The third source of variability in the dynamic performance monitoring related to the specific target IC within a subject. Across different monitored ICs within the same subject, the results of Table 1 and Figure 2 confirm that the difference in behavior of different subjects was consistent across ICs. Indeed the performance improved for all the subjects when monitoring task related RVMT and LVMT (Figure 1) with respect to Resting State Network (RSN) related Default Mode Network (DMN) (Figure 2) thus proving that difference in the nature of monitored IC was the strongest source of variability (up to 30%) for these kind of presented methods. This may be due to the fact that different activations have particular statistical distribution properties, being more or less suitable to be extracted by ICA algorithms. In addition, a reason for the difficulty in extracting spatial characterization of the DMN is its low frequency relative to the sliding time window length, thus making it difficult for the algorithm to correctly follow it.

This work has some limitations. One limitation is related to the definition of dynamic monitoring performance, which depends on temporal or spatial correlations with a template reference derived from the whole time course. It is not necessarily correct to expect that spatial-temporal characteristics derived from the sliding window along the time course should match the ones derived from the whole time course. For this reason the performance measures are only indicative.

The possibility that the actual dynamic brain activation is correctly identified by these on-line methods opens the door to future definition of techniques and experiments. These experiments could exploit these methods to have a confirmation of transient activation identification independently of the off-line analysis, which represents a general reference for evaluation of results, but may also represent a bias.

Another limitation relates to the simulation nature of the work, which should be further evaluated on a real implementation in which the performance of the different methods can be studied, for example using a neurofeedback setup.

The comparison of ICA with non-ICA approaches in rt-fMRI setups was beyond the scope of this study. Given the known potential advantages of data-driven analysis (Norman et al., 2006; Magland et al., 2011; Sitaram et al., 2011), and in particular ICA methods (Esposito et al., 2003; Beckmann and Smith, 2004; Calhoun et al., 2009), this simulation study is limited to the comparison of novel ICA-based methods for rt-fMRI using robust activation in well-known visual motor areas. Future studies will be needed to evaluate these ICA methods with brain activation that could be more challenging to identify.

This work proposed and evaluated several strategies for using a priori information for the monitoring of brain networks in real-time fMRI experiments. The performance of the methods was characterized by both computation speed and correlation between the spatial-temporal properties of a target independent component derived dynamically and a reference component. The method that gave the highest performance was based on the Back-projection of a constant target spatial map derived by the spatial localizer. In this method the use of ICA was exploited only in the Functional Localizer phase, while during the on-line monitoring the reference component was kept constant and not updated with any ICA algorithm. This combination of both ICA and non-ICA methods shows thus to be very helpful and promising. This method had the limitation that its reference was constant and this means that it may not be optimal to follow dynamic changes as it cannot adapt to changes in brain. The other tested methods were based on the use of adaptive spatial, temporal, or spatial-temporal priors and may have useful applications in studies where there is a need of higher flexibility to monitor variable activation.

Conclusion

This study proposed and evaluated several strategies for using spatial and/or temporal a priori information in ICA-based methods for the monitoring of brain networks in real-time fMRI experiments. The effectiveness of the novel real-time ICA-based method was evaluated against the off-line ICA analysis that is typically possible after all data has been acquired. The performance of the methods was characterized by both computation speed and correlation between the spatial-temporal properties of a target independent component derived dynamically and a reference component. In our testing conditions of relative low frequency task-induced activations (with a period of 20 s) we found that the Back-projection method outperformed the other methods giving the highest spatial-temporal correlations to the reference and the fastest computation time. The Back-projection method here investigated uses the ICA decomposition only in the functional localizer data, and not during the dynamic on-line analysis. It remains to be further investigated whether the spatial-temporal constrained methods can be better, as in principle expected, in situations where the networks to be monitored have higher frequency fluctuations in space and time.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work has been realized thanks to the support from the Provincia autonoma di Trento and the Fondazione Cassa di Risparmio di Trento e Rovereto, Italy. Jorge Jovicich was supported by a PRIN grant (Prot. 2008HETS3W _002) from MIUR (Italy). The authors would also like to specially thank James Hartzell.

References

Anderson, A., Han, D., Douglas, P. K., Bramen, J., and Cohen, M. S. (2011). “Real-time functional MRI classification of brain states using markov-svm hybrid models: peering inside the rt-fMRI black box,” in MLINI. Machine Learning and Interpretation in Neuroimaging – International Workshop, MLINI 2011, Held at NIPS 2011. Lecture Notes in Computer Science, Vol. 7263, eds G. Langs, I. Rish, M. Grosse-Wentrup, and B. Murphy (Sierra Nevada: Springer), 242–255. [December 16–17, 2011; Revised Selected and Invited Contributions].

Beckmann, C. F., and Smith, S. M. (2004). Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23, 137–152.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bell, A. J., and Sejnowski, T. J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Comput. 7, 1129–1159.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calhoun, V., and Adali, T. (2006). Unmixing fMRI with independent component analysis. IEEE Eng. Med. Biol. Mag. 25, 79–90.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calhoun, V. D., Adali, T., Pearlson, G. D., and Pekar, J. J. (2001). Spatial and temporal independent component analysis of functional MRI data containing a pair of task-related waveforms. Hum. Brain Mapp. 13, 43–53.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calhoun, V. D., Eichele, T., and Pearlson, G. (2009). Functional brain networks in schizophrenia: a review. Front. Hum. Neurosci. 3:17. doi:10.3389/neuro.09.017.2009

CrossRef Full Text

Caria, A., Sitaram, R., and Birbaumer, N. (2012). Real-time fMRI: a tool for local brain regulation. Neuroscientist 18, 487–501.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chiu, P., Lisinski, J., King-Casas, B., Sharp, J., Eagleman, D., Versace, F., et al. (2010). Modulating “crave” and “don’t crave” brain states with real-time fMRI neurofeedback in smokers. Paper Presented at the 16th Annual Meeting of the Organization for Human Brain Mapping, Barcelona.

Cox, R. W., Jesmanowicz, A., and Hyde, J. S. (1995). Real-time functional magnetic resonance imaging. Magn. Reson. Med. 33, 230–236.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

deCharms, R. C. (2008). Applications of real-time fMRI. Nat. Rev. Neurosci. 9, 720–729.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

deCharms, R. C., Maeda, F., Glover, G. H., Ludlow, D., Pauly, J. M., Soneji, D., et al. (2005). Control over brain activation and pain learned by using real-time functional MRI. Proc. Natl. Acad. Sci. U.S.A. 102, 18626–18631.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

DeMartino, F., Gentile, F., Esposito, F., Balsi, M., Salle, F. D., Goebel, R., et al. (2007). Classification of fMRI independent components using IC-fingerprints and support vector machine classifiers. Neuroimage 34, 177–194.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Esposito, F., Seifritz, E., Formisano, E., Morrone, R., Scarabino, T., Tedeschi, G., et al. (2003). Real-time independent component analysis of fMRI time-series. Neuroimage 20, 2209–2224.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gembris, D., Taylor, J. G., Schor, S., Frings, W., Suter, D., and Posse, S. (2000). Functional magnetic resonance imaging in real time (fire): sliding-window correlation analysis and reference-vector optimization. Magn. Reson. Med. 43, 259–268.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Goebel, R. (2012). Brainvoyager – past, present, future. Neuroimage 62, 748–756.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Himberg, J., Hyvrinen, A., and Esposito, F. (2004). Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage 22, 1214–1222.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hinds, O., Ghosh, S., Thompson, T. W., Yoo, J. J., Whitfield-Gabrieli, S., Triantafyllou, C., et al. (2011). Computing moment-to-moment bold activation for real-time neurofeedback. Neuroimage 54, 361–368.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hyvrinen, A., and Oja, E. (2000). Independent component analysis: algorithms and applications. Neural Netw. 13, 411–430.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

LaConte, S. M. (2011). Decoding fMRI brain states in real-time. Neuroimage 56, 440–454.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

LaConte, S. M., Peltier, S. J., and Hu, X. P. (2007). Real-time fMRI using brain-state classification. Hum. Brain Mapp. 28, 1033–1044.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lin, Q. H., Liu, J., Zheng, Y. R., Liang, H., and Calhoun, V. D. (2010). Semiblind spatial ICA of fMRI using spatial constraints. Hum. Brain Mapp. 31, 1076–1088.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Maclaren, J., Herbst, M., Speck, O., and Zaitsev, M. (2013). Prospective motion correction in brain imaging: a review. Mag. Reson. Med. 69, 621–636.

CrossRef Full Text

Magland, J. F., Tjoa, C. W., and Childress, A. R. (2011). Spatio-temporal activity in real time (star): optimization of regional fMRI feedback. Neuroimage 55, 1044–1053.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

MATLAB. (2010). Version 7.10.0 (R2010a). Natick, MA: The MathWorks Inc.

McKeown, M. J., Makeig, S., Brown, G. G., Jung, T. P., Kindermann, S. S., Bell, A. J., et al. (1998). Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6, 160–188.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mouro-Miranda, J., Bokde, A. L. W., Born, C., Hampel, H., and Stetter, M. (2005). Classifying brain states and determining the discriminating activation patterns: support vector machine on functional MRI data. Neuroimage 28, 980–995.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Norman, K. A., Polyn, S. M., Detre, G. J., and Haxby, J. V. (2006). Beyond mind-reading: multi-voxel pattern analysis of fMRI data. Trends Cogn. Sci. (Regul. Ed.) 10, 424–430.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Posse, S., Fitzgerald, D., Gao, K., Habel, U., Rosenberg, D., Moore, G. J., et al. (2003). Real-time fMRI of temporolimbic regions detects amygdala activation during single-trial self-induced sadness. Neuroimage 18, 760–768.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shibata, K., Watanabe, T., Sasaki, Y., and Kawato, M. (2011). Perceptual learning incepted by decoded fMRI neurofeedback without stimulus presentation. Science 334, 1413–1415.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sitaram, R., Lee, S., Ruiz, S., Rana, M., Veit, R., and Birbaumer, N. (2011). Real-time support vector classification and feedback of multiple emotional brain states. Neuroimage 56, 753–765.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Soldati, N., Calhoun, V. D., Bruzzone, L., and Jovicich, J. (2013). ICA analysis of fMRI with real-time constraints: an evaluation of fast detection performance as function of algorithms, parameters and a priori conditions. Front. Hum. Neurosci. 7:19. doi:10.3389/fnhum.2013.00019

CrossRef Full Text

Thompson, T., Hinds, O., Ghosh, S., Lala, N., and Triantafyllou, C. (2009). Training selective auditory attention with real-time fMRI feedback. Neuroimage 47, S65.

CrossRef Full Text

Weiskopf, N. (2012). Real-time fMRI and its application to neurofeedback. NeuroImage 62, 682–692.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weiskopf, N., Scharnowski, F., Veit, R., Goebel, R., Birbaumer, N., and Mathiak, K. (2004). Self-regulation of local brain activity using real-time functional magnetic resonance imaging (fMRI). J. Physiol. Paris 98, 357–373.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weiskopf, N., Sitaram, R., Josephs, O., Veit, R., Scharnowski, F., Goebel, R., et al. (2007). Real-time functional magnetic resonance imaging: methods and applications. Magn. Reson. Imaging 25, 989–1003.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text