Subject order-independent group ICA (SOI-GICA) for functional MRI data analysis

NeuroImage 51 (2010) 1414–1424

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g

Subject order-independent group ICA (SOI-GICA) for functional MRI data analysis Han Zhang a, Xi-Nian Zuo b, Shuang-Ye Ma a, Yu-Feng Zang a, Michael P. Milham b, Chao-Zhe Zhu a,⁎ a b

State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, Beijing 100875, China Phyllis Green and Randolph Cōwen Institute for Pediatric Neuroscience at the New York University Child Study Center, New York, USA

a r t i c l e

i n f o

Article history: Received 11 February 2010 Accepted 12 March 2010 Available online 23 March 2010 Keywords: Independent component analysis (ICA) Group inference Inconsistency Variability Instability Functional magnetic resonance imaging (fMRI) Resting state

a b s t r a c t Independent component analysis (ICA) is a data-driven approach to study functional magnetic resonance imaging (fMRI) data. Particularly, for group analysis on multiple subjects, temporally concatenation group ICA (TC-GICA) is intensively used. However, due to the usually limited computational capability, data reduction with principal component analysis (PCA: a standard preprocessing step of ICA decomposition) is difficult to achieve for a large dataset. To overcome this, TC-GICA employs multiple-stage PCA data reduction. Such multiple-stage PCA data reduction, however, leads to variable outputs due to different subject concatenation orders. Consequently, the ICA algorithm uses the variable multiple-stage PCA outputs and generates variable decompositions. In this study, a rigorous theoretical analysis was conducted to prove the existence of such variability. Simulated and real fMRI experiments were used to demonstrate the subjectorder-induced variability of TC-GICA results using multiple PCA data reductions. To solve this problem, we propose a new subject order-independent group ICA (SOI-GICA). Both simulated and real fMRI data experiments demonstrated the high robustness and accuracy of the SOI-GICA results compared to those of traditional TC-GICA. Accordingly, we recommend SOI-GICA for group ICA-based fMRI studies, especially those with large data sets. © 2010 Elsevier Inc. All rights reserved.

Introduction Independent component analysis (ICA) was originally proposed as a blind source separation algorithm in signal processing and has proven a powerful, efficient and reliable tool for exploratory analyses in functional magnetic resonance imaging (fMRI) studies (McKeown et al., 1998; Kiviniemi et al., 2003; Hu et al., 2005; Zuo et al., 2010). Early applications of ICA focus on single-subject fMRI data (i.e., individual-level) (McKeown et al., 1998). To address the challenge of matching single-subject components between subjects in group-level analyses, group ICA (GICA) was developed (see Guo and Pagnoni, 2008 and Calhoun et al., 2009 for reviews). Temporally concatenation GICA (TC-GICA) has emerged as one of the leading GICA approaches (Calhoun et al., 2001; Beckmann et al., 2005). TC-GICA concatenates all individual data along the temporal dimension before ICA decomposition. However, after fMRI data are concatenated across subjects, the size of the merged data usually become quite large. Consequently, as a standard preprocessing step of ICA, principal component analysis (PCA) is difficult to achieve due to limited computational capability, especially for large group studies. To solve this problem, a hierarchically multi-stage PCA (MS-PCA) data reduction is

⁎ Corresponding author. Fax: +86 10 58806154. E-mail address: [email protected] (C.-Z. Zhu). 1053-8119/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2010.03.039

often adopted and has been implemented in GIFT (Calhoun et al., 2001) and MELODIC (Smith et al., 2004). Briefly, after the first stage individuallevel PCA reduction, additional stages of PCA reductions are employed to further reduce the data dimension in a hierarchical manner. In general, the larger the concatenated data are, the more PCA-reduction stages are required. Currently, TC-GICA with three-stage MS-PCA reduction (TC-GICA3) is widely used, and several interesting findings have been reported (Calhoun et al., 2008; Jafri et al., 2008; Chen et al., 2008; Stevens et al., 2009). However, TC-GICA3 has a potentially hazardous caveat, as it may produce different results due to different concatenation orders of the individual fMRI data sets. Accordingly, much attention should be paid to the conclusions drawn based on results from TC-GICA3. Indeed, it is necessary to systematically investigate such subject concatenation order (SCO)-related variability of results and to propose an algorithm to overcome this drawback. In this study, we first prove the existence of SCO-induced variability in three-stage MS-PCA reduction results and in turn the following ICA decompositions (i.e., TC-GICA3 results). Furthermore, well-designed simulations were used to systemically demonstrate how various parameters in TC-GICA3 affect the SCO-induced variability. Finally, a subject order-independent group ICA (SOI-GICA) was proposed to solve this problem. Both simulated and real fMRI experiments were conducted to evaluate the performance of the proposed SOI-GICA algorithm, and some limitations of the proposed SOI-GICA algorithm are also discussed.

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

Theoretical analyses

Formulations of MS-PCA3 data reduction Here, we provide formulized forms of MS-PCA3 (see Appendix A for the formulations of MS-PCA1 and MS-PCA2). Suppose we have N subjects. The first stage of PCA reduction is conducted at the individual level. The second stage of PCA reduction utilizes a subjectgrouping strategy, where every S (S b N) reduced individual data sets are successively grouped, and the whole N data sets are equally divided into G subgroups. Note that the equal division of total N subjects simplifies notation, but the ideas can be easily generalized to unequal subgroup size. The concatenated data for each subgroup are reduced at the subgroup level: 3 ð1;g;1Þ X 6 MSPCA3 7 6 ð1;g;2Þ 7 6 XMSPCA3 7 7 6 7
−1 6 ⋮ ð2;gÞ 7 6 = FMSPCA3 7; 6 6 X ð1;g;sÞ 7 6 MSPCA3 7 7 6 7 6 ⋮ 5 4 ð1;g;SÞ XMSPCA3 2

ð2;gÞ

S CNS CN−S :::CSS different G PG grouping manners. Different grouping manners result in different inputs to the second-stage PCA reduction and different data reduction outputs. The different outputs, in turn, lead to different MS-PCA3 outputs, X(3) MS − PCA3. A rigorous mathematical proof on “different grouping manners produce different PCA results” is given in Appendix B. Of note, different SCOs may belong to the same grouping manner. For example, if we switch the subject-concatenation order within a subgroup, all possible SCOs yield only one grouping manner. We named the SCOs that will yield distinct grouping manners as “effective” SCOs. C S C S :::CSS effective SCOs in total N! different SCOs. Thus, there are N N−S PGG

According to combinatorial theory, there are

For simplicity, we refer to TC-GICA with an n-stage MS-PCA reduction (MS-PCAn) as TC-GICAn (e.g., TC-GICA1 denotes TC-GICA with MS-PCA1, and TC-GICA3 is TC-GICA with MS-PCA3).

XMSPCA3

1415

ð2;1Þ

XMSPCA3

ð3Þ

ð1Þ

3

7 6 7
−1 6 ð2;2Þ ð3Þ ð3Þ 6 XMSPCA3 7 XMSPCA3 = FMSPCA3 7: 6 7 6 ⋮ 5 4

(named as “PC-TCTC-GICA3 uses MS-PCA3 outputs, X(3) MS − PCA3 GICA3”), as the input of ICA decomposition: XMSPCA3 = Aˆ TCGICA3 SˆTCGICA3 ;

(1,g,s) where XMS denotes the reduced data from the first-stage of PCA − PCA3 reduction for subject s in subgroup g (g = 1, 2,..., G). The third stage of (2,g) (g = 1, PCA reduction is applied to the concatenated data of XMS − PCA3 2,..., G):

2

TC-GICA3 outputs depending on SCO

ð2Þ

ð2;GÞ

XMSPCA3

ð3Þ

where ŜTC − GICA3 denotes the estimated spatial maps from "aggregated" data (named as “aggIC-TC-GICA3”). It is quite straightforward that may cause variable ICA decomposition. Consequentvariable X(3) MS − PCA3 ly, the reconstructed (Calhoun et al., 2001; Beckmann et al., 2009) individual-level components (“indivIC-TC-GICA3”) and the following group-level statistics (“groupStats-TC-GICA3”) may also change. In brief, the dependence of ICA decomposition on SCO originates from the dependence of MS-PCA3 reduction on SCO. Fig. 1 depicts the origination and propagation of the SCO-dependence in a whole TC-GICA3 algorithm. From the viewpoint of matrix eigenvalue decompositions, MSPCA1 uses covariance information from all variables (the whole concatenated data). However, MS-PCA3 uses covariance information from only parts of all variables in a hierarchical manner. This makes MS-PCA3 an approximation of MS-PCA1. Thus, TC-GICA3 results are approximations of TC-GICA1 results. The approximated error is defined as one minus the absolute value of Pearson's r between the TC-GICA3 results (θTC − GICA3) and the TC-GICA1 results (θTC − GICA1):

MS-PCA3 output may change with different SCOs

εθTCGICA3 = 1−jρθTCGICA3 ;θTCGICA1 j:

When conducting MS-PCA3 reduction, the grouping manner in second-stage PCA reduction may change with different SCO. Specifically, consider N subjects to be equally divided into G subgroups with S subjects in each subgroup in the second stage of PCA reduction.

Of note, εθTC-GICA3 depends on SCO. The variability in εθTC-GICA3 for multiple SCOs is thus defined as a measurement of the SCOdependence of TC-GICA3 outputs.

ð4Þ

Fig. 1. Origin and propagation of SCO-induced variability. Different subgrouping manner in subgroup formation in the second stage PCA reduction introduces SCO-related variability (indicated by a red arrow) in the reduced data. Such variability propagates through the following ICA decomposition, back-reconstruction (BR), and group analysis (GA), resulting in variable aggregate component (Aggregate IC), individual-level component (Individual IC), and group-level statistics (Group stats). The figure also shows randomized IV-induced variability (indicated by another red arrow).

1416

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

SOI-GICA algorithm To address the problem of SCO-dependence in TC-GICA3, “Subject Order-Independent Group ICA” (SOI-GICA) is proposed and consists of three steps: Step 1 Generating multiple different effective SCOs. We generated effective SCOs by randomly permuting the natural order of subjects and removing the redundant ones; Step 2 Performing TC-GICA3 with each effective SCO; and Step 3 Integrating multiple TC-GICA3 outputs. According to the Law of Large Numbers (Grimmett and Stirzaker, 1992), when integrating multiple TC-GICA3 results, SOI-GICA may improve the approximation of the TC-GICA1 result and thus lead to a SCO-independent result. Specifically, with T different effective SCOs, TC-GICA3 generates T results. SOI-GICA integrates all of the T results and produces the corresponding results: i. AggIC-SOI-GICA is calculated by averaging multiple aggIC-TCGICA3 results according to Eq. (5). Of note, this result is provided as an intermediate result: 1 T SˆSOIGICA = ∑ SˆTCGICA3 ðtÞ T t =1

ð5Þ

n ii. IndivIC-SOI-GICA (spatial maps ŜSOI − GICA and the associated n ) for subject n (n = 1, 2,..., N) are calculated ̂ time courses T SOI-GICA by averaging across multiple indivIC-TC-GICA3 results:

f

g

1 T n n Sˆ SOIGICA = ∑ SˆTCGICA3 ðtÞ T t=1 1 T n n Tˆ SOIGICA = ∑ Tˆ TCGICA3 ðtÞ T t =1

Experiments Aside from the theoretical analysis, three simulated experiments were designed to demonstrate the existence of SCO-dependence in TC-GICA3 and the improved performance of SOI-GICA. First, to demonstrate the existence of SCO-dependence, a multi-dimension scaling (MDS) approach (Mead, 1992) was utilized to project the outputs from multiple TC-GICA3 with different SCOs in a highdimensional space onto a two-dimensional plane. Then, to quantitatively evaluate the improved performance of SOI-GICA, εθSOI-GICA and εθTC-GICA3 were compared. In addition, we investigated the influence of various parameters mentioned in the SOI-GICA algorithm section on the performance of SOI-GICA via stepwise varying one parameter while keeping others fixed. A general flowchart of the simulated experiments is depicted in Fig. 2.

Simulated experiment 1: task-condition within-group analysis ð6Þ

iii. GroupStats-SOI-GICA is generated by statistical analysis of indivIC-SOI-GICA across all subjects. Due to the fact that SOI-GICA is based on multiple TC-GICA runs, the underlying assumption of the SOI-GICA is that all subjects share the same spatial maps without the requirement of the same time courses (i.e., the same assumption as TC-GICA (Calhoun et al., 2001)). Similar to Eq. (4), the approximation error of the SOI-GICA (εθSOI-GICA) to the TC-GICA1 is defined as: εθSOIGICA = 1−jρθSOIGICA ;θTCGICA1 j;

i. the parameters in subgroup formation and the following second-stage PCA reduction, i.e., the subgroup size (S) and the retained dimension after the second-stage PCA reduction (J); ii. the parameter in ICA estimation, i.e., the initial value (IV); iii. the parameters related to dataset itself, i.e., the contrast-tonoise rate (CNR); and iv. the convergence of the SOI-GICA algorithm, i.e., the number of TC-GICA3 runs required for SOI-GICA (T).

ð7Þ

Additionally, we systemically investigated the possible impact of the following parameters on the performance of SOI-GICA:

Data sets for a group of 30 subjects were generated to simulate task-condition fMRI data for within-group analysis (TG1), according to the generation model:

A

A

B

B

Di = T × S + T × S + Ni :

ð8Þ

In Eq. (8), Di represents the data set for subject i (i = 1, 2,..., 30), SA and SB represent two spatial maps related to task A and B, and TA and TB denote the corresponding task-evoked temporal activities (TA is three times stronger than TB). Both T and S are identical for all subjects (as shown in Fig. 3(a)). The sign “×” denotes the exterior product, and Ni represents the additive subject-specific Gaussian noise (CNR = 3.9). For additional details, please see the Supplementary materials. Using both conventional TC-GICA3 and our proposed SOI-GICA, the individual- and group-level spatial maps related to both task A and B were estimated and compared.

Fig. 2. Flowchart of the simulated experiments. This flowchart summarizes the main steps of the simulated experiments in this paper, including the calculations of TC-GICA3 (indicated by a dotted rectangle), SOI-GICA, and TC-GICA1 (bottom, as a referential baseline). SOI-GICA integrates multiple TC-GICA3 runs with different effective SCOs (from 1 to T).

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

1417

Fig. 3. Sources for three simulated experiments. (a) Two sources for the task-condition within-group analysis (TG1); (b) three sources of a sample subject i for the resting-state within-group analysis (RG1); and (c) two sources of two groups (Group 1 and Group 2) for the task-condition between-group analysis (TG2). S represents the spatial map and T the time course.

Simulated experiment 2: resting-state within-group analysis

Simulated experiment 3: task-condition between-group analysis

Resting-state data sets for a group of 30 subjects (RG1) were simulated according to the model:

To investigate the SCO-dependence in group difference analysis, which is of interest in various psychiatric and neurological disorder studies (Greicius et al., 2004), task-condition data sets from two groups (TG2) were generated according to:

A

A

B

B

C

C

Di = Ti × S + Ti × S + Ti × S + Ni :

ð9Þ

In Eq. (9), the spatial map SA is designed to simulate a functional system of interest, with subject-specific low-frequency spontaneous fluctuation TiA (i = 1, 2,..., 30). The spatial maps S B and S C represent two noise-related components with associated subject-specific time courses TiB and TiC. Fig. 3(b) shows the three spatial maps and the associated time courses for a subject (for more details, please see the Supplementary materials). As low-frequency spontaneous fluctuation is usually a research interest in resting-state fMRI studies, the individual-level time course TiA and group-level spatial map S A were used to demonstrate results.

A

A

B

B

Dij = T × S + T × Sj + Nij :

ð10Þ

In Eq. (10), Dij and Nij represent the simulated data and noise (CNR = 1) for subject i (i = 1, 2,..., 15) in group j (j = 1, 2). Group difference was simulated to only lie in SB (Fig. 3(c)) (for more details, please sees the Supplementary materials). The group difference in SB was assessed by a two-sample t-test on the individual-level spatial maps. The derived group difference (t-maps) from TC-GICA3 and SOIGICA were compared.

1418

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

Real resting-state fMRI experiment A group of resting-state BOLD fMRI data sets from 20 healthy subjects (10 males, 22.80 ± 3.65 years old) were acquired by a Siemens 3.0T Trio scanner at Beijing Normal University. The data preprocessing was conducted using SPM5 (http://www.fil.ion.ucl.ac. uk/spm/), including acquisition-time correction, head-motion correction, spatial normalization, and smoothing (the detailed information can be found in the Supplementary materials). In each TC-GICA3 run, the component number was estimated to be 53 using a maximum description length approach (Calhoun et al., 2001; Li et al., 2007). MS-PCA3 was performed with 53 principal components retained for all reduction stages. ICA decomposition was conducted using FastICA (Hyvarinen, 1999) with the following options: randomized initial value, symmetric approach, nonlinearity= pow3, maximum number of iterations= 1000, fine-tune off, epsilon = 0.0001, and stabilization = off. The individual-level components were obtained from back-reconstruction and converted into z-scores. The group-level results were generated by one-sample t-test. The SOI-GICA result was derived based on the outputs from 100 TC-GICA3 runs with different effective SCOs, according to convergence studies in the simulated experiments (see Results). Results Simulated experiment 1: task-condition within-group analysis TC-GICA3 was carried out 100 times with different effective SCOs (the other parameters were: S = 5, J = 10, and fixed initial value in FastICA algorithm). The results from each processing stage, including

outputs from MS-PCA3, aggIC-, inidivIC-, and groupStats-TC-GICA3, are shown in Fig. 4(a)–(e), respectively, using MDS plotting. Each dot represents the result from a single TC-GICA3 run. The scattered patterns in the two-dimensional plane clearly demonstrates SCOinduced variability. Moreover, such variability arose after MS-PCA3 reduction and sequentially propagated to the following ICA decomposition, individual components reconstruction, and the final group statistics. Of note, as another potential source of ICA instability (Ylipaavalniemi and Vigario, 2008), the randomized initial value (IV) was controlled by fixing it to the same value (i.e., the identity matrix) for all 100 TC-GICA3 runs. The impact of S (i.e., the subgroup size in second-stage PCA) on the performance of TC-GICA3 (dots) and SOI-ICA (circles) is shown in Fig. 5(a) and (b) for estimated SA at both individual- (left) and group- (right) level. For each S value, both the indivIC- and groupStats-SOI-GICA exhibited better approximation than nearly all of the corresponding TC-GICA3 results. Moreover, the mean and variance of εθTC-GICA3 across 100 runs were enlarged as S decreased from 15 to 2, indicating greater SCO-dependence for smaller subgroup size. Conversely, SOI-GICA generated quite stable results, indicating that it is largely unaffected by the parameter of subgroup size. For the other source, SB, the SCO-dependence of the indivIC- and groupStats-TCGICA3 was more serious due to its weaker activation (CNR = 1.3) than that for SA(CNR = 3.9) (Supplementary Fig. 1). Taking indivIC-TCGICA3 for example, the approximated errors of SB across 100 runs varied from ∼0.15 to ∼0.37 (i.e., the Pearson's r ≈ 0.63–0.85), while those for SA were only from ∼0.03 to ∼0.09. SOI-GICA, however, yielded much more robust and accurate results at both individualand group-level, even in the situation of weak activation. It should be noted that we mainly presented the results of component SA because

Fig. 4. Illustration of the TC-GICA3 variability induced by different SCOs. TC-GICA3 was performed 100 times on simulated data TG1, each time with a different effective SCO but a constant IV. Each cross denotes the estimated source A from a single TC-GICA3 run. Five MDS plots represent the variable results at all processing stages of TC-GICA3, including MS-PCA3 (a), aggregate IC (b), spatial map and time course of the individual-level component (c, d) and group-level statistics (e).

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

1419

Fig. 5. Comparison of accuracy and robustness between SOI-GICA (circles) and TC-GICA3 (dots) (TG1, source A). Left panels: individual-level spatial map for a sample subject; Right panels: group-level statistical map. First row: S (i.e., subgroup size) effect; second row: J (i.e., number of retained dimensions after second-stage PCA reduction) effect.

the results of component SB were highly similar. Thus, in the remainder of the article, we also present only part of the results for the same reason. The impact of J (i.e. the number of retained dimensions after second-stage PCA reduction) on the performance of TC-GICA3 and SOI-GICA is shown in Fig. 5(c) and (d) in a similar way as that for parameter S. SOI-GICA clearly yielded more accurate individual- and group-level results than TC-GICA3 at each J (from 2 to 90). Moreover, SOI-GICA generated more robust results to J than TC-GICA3. Randomization of the initial value (IV) in ICA algorithms (e.g., FastICA) induces variability in ICA decomposition (Ylipaavalniemi and Vigario, 2008). Considering this, questions are raised as to which of the two variability sources (i.e., IV and SCO) is dominant in TC-GICA3 and whether the proposed SOI-GICA still works well with the additional variability induced by randomized IV. Fig. 6 shows the variability independently contributed by different effective SCOs (black dots) and that by randomized IV (blue dots). The difference between the maximal and the minimal εθTC-GICA3 across 100 TC-GICA3 runs was adopted to quantitatively evaluate such variability. For the individual-level spatial map of SA(S = 2 and J = 10), the randomized IV-induced variability was 0.0001, while the different effective SCOinduced variability was 0.0383 (∼ 340 times more). Fig. 6 also shows

the variability contributed simultaneously by IV and SCO (red dots). Although one additional source of variability (i.e., IV) was introduced, the variability (0.0387) was almost the same as that contributed only by the SCO (0.0383). This result indicated that SCO is a dominant source of variability in TC-GICA3. Moreover, SOI-GICA produced results with almost the same accuracy, even when the additional variability source (IV) was present (see black and red circles). The robustness of the SOI-GICA algorithm was evaluated under different data CNR (from 5.5 to 2). As data CNR decreased, both εθTC-GICA3 and εθSOI-GICA for spatial maps of SA were enlarged. However, SOI-GICA produced more accurate results than most of the TC-GICA3 results at all CNR levels (Supplementary Fig. 2(a) and (b)). The basic idea behind SOI-GICA is to integrate the outputs from multiple TC-GICA3 with different SCOs to generate a more accurate result. It is thus important to investigate the convergence property of the SOI-GICA algorithm and to determine a proper number of different effective SCOs (i.e., T). As T increased from 2 to 2000, integrating more TC-GICA3 results, SOI-GICA tended to yield increasingly accurate results (Fig. 7). This is consistent with the Law of Large Numbers. Specifically, as T reached ∼ 50, the individual-level SOI-GICA result converged (Fig. 7(a)). For the group-level result (Fig. 7(b)), SOIGICA converged more quickly (T = ∼10).

1420

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

Simulated experiment 3: task-condition between-group analysis As shown in Fig. 8(a), the approximated errors of the groupdifference detected by 100 TC-GICA3 runs with different effective SCOs and fixed IV were dispersed over a large range (≈ 0.3 - 0.9). This indicated that while several TC-GICA3 results were acceptable, several others were inaccurate. In contrast, SOI-GICA outperformed TC-GICA3 with a much lower approximated error of 0.2. The peak t-value in group-difference map as shown in Fig. 8(b) also demonstrated the improved performance (i.e., the higher t values) of SOI-GICA. Moreover, the ROC curves in Fig. 8(c) further confirmed that SOIGICA detected group differences with much better sensitivity and specificity (red). In contrast, TC-GICA introduced more false-positive and false-negative detections (black). The convergence property of SOI-GICA for group-difference detection was investigated. With each S setting from 2 to 15, the group comparison result from SOI-GICA always became stable before T reached 100 (Supplementary Fig. 5). Fig. 6. Comparison between variability induced by two sources: SCO and IV. First column: 100 TC-GICA3 results (dots) and the corresponding SOI-GICA result (circle) with different effective SCOs but fixed IV; second column: those results with randomized IV but fixed SCO; and third column: those results with different effective SCOs and randomized IVs.

Simulated experiment 2: resting-state within-group analysis Considering the limited space, we mainly focus on subject-specific individual-level time course of TAi in this “resting-state” experiment. Also, we only demonstrate the impacts of parameters S and J on the performance of SOI-GICA and TC-GICA3. We fixed IV across multiple TC-GICA3 runs in order to rule out the potential effect of randomized IV (although its effect can be ignored as demonstrated in Simulated experiment 1: task-condition within-group analysis). As shown in Supplementary Fig. 3(a) and (c), the individual-level time course TAi derived from SOI-GICA was more accurate than TCGICA3 at nearly every S (from 2 to 15) and J (from 4 to 80) value. Moreover, the accuracy of estimated TAi from SOI-GICA was robust to these parameters. The improved performance was also demonstrated by SOI-GICA-derived group-level spatial map SA(Supplementary Fig. 3 (b) and (d)). The convergence property of the SOI-GICA algorithm was also investigated. Our results demonstrated that the SOI-GICA results became stable as T approached ∼100 (Supplementary Fig. 4).

Real resting-state fMRI experiment In a real fMRI experiment, SOI-GICA was carried out with randomized IV to demonstrate its performance in a more practical sense. Five intensively studied resting-state networks (RSNs) involving both primary and higher functioning were chosen to demonstrate the results by using a template-matching approach in GIFT. The prior spatial templates were generated based on another independent resting-state fMRI data set (Lv et al., 2008). These templates were networks associated with visual (VIS), sensorimotor (MOT), and auditory processing (AUD), as well as the default-mode network (DMN) and executive control network (EXE). Consistent with findings from previous resting-state fMRI studies (Beckmann et al., 2005; Damoiseaux et al., 2006, 2008; Kiviniemi et al., 2009; Smith et al., 2009; Zuo et al., 2010), the groupStats-SOIGICA of the five components of interest are shown in the first row of Fig. 9. For comparison, the two instances of 100 groupStats-TC-GICA3 results are shown in the second and third rows, respectively. When compared with SOI-GICA, TC-GICA3 was likely to either fail to detect the brain areas which should belong to those RSNs or to detect irrelevant brain areas by mistake (as indicated by arrows). Taking the MOT component as an example, only 41% of the TC-GICA3 runs detected the nucleus in the mesencephalon bilaterally; while 27%

Fig. 7. SOI-GICA convergence (TG1, source A). The approximated error of the SOI-GICA result (left: individual-level spatial map for a sample subject; right: group-level statistical map) was plotted against the number of integrated TC-GICA3 runs with different effective SCOs (i.e., T).

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

1421

Fig. 8. Comparisons of group-difference detection performance between SOI-GICA (circles) and TC-GICA3 (dots). (a) The result accuracy (measured by approximated error to the referential baseline); (b) the peak t-value in group-difference map; and (c) the ROC curves for the two approaches (SOI-GICA in red and TC-GICA3 in black).

of the runs made one-sided detection, and 32% failed in detection (see more detailed comparisons in Supplementary materials). The relationship between SOI-GICA and TC-GICA3 results is depicted in MDS plots both for individual- (Supplementary Figs. 6–10) and grouplevel (Fig. 9) results. All MDS plots demonstrate a notable variability of TC-GICA3 results. The degree of such variability (indicated by the scattering range of all TC-GICA3 results in the plot) was generally larger for RSNs involved in higher cognitive functioning than those involved in primary functioning.

Discussion In this study, SCO-induced TC-GICA3 variability was demonstrated using the theoretical analysis, as well as the simulated and real fMRI experiments. Specifically, different grouping manners in subgroup formation during the second-stage PCA reduction introduces SCOrelated variability in MS-PCA3 outputs. Such variability propagates to ICA decomposition and subsequent processing that uses MS-PCA3 reduced data as an input. In particular, we noted that the SCO-induced variability in group-difference analysis was more serious (increased by 25-fold, given the same noise level) than that in within-group analysis. Those results suggest that the findings from TC-GICA3 with a specified SCO should be carefully interpreted. We speculated that such variability might, in all likelihood, be present in following analyses of TC-GICA3 outputs, such as correlation analysis between TC-GICA3 results and behavior or clinical measures (Damoiseaux et al., 2008), reliability analysis (Chen et al., 2008), and “functional network connectivity” analysis (Jafri et al., 2008). Aside from SCO-induced variability, randomized IV also introduces variability in TC-GICA3. Indeed, this variability source has been widely investigated in individual-level ICA studies (Duann et al., 2003; Himberg et al., 2004; Yang et al., 2008; Ylipaavalniemi and Vigario, 2008). The way to deal with IV-induced variability is to average results from multiple ICA decompositions with randomized IVs to generate IV-independent outputs (Himberg et al., 2004, Yang et al., 2008; Ylipaavalniemi and Vigario, 2008). Through numerical simulations, we concluded that the IV-induced variability is relatively much smaller than the SCO-induced variability. Thus, in TC-GICA3, SCO is the dominant variability source and only investigating IV-induced variability in TC-GICA3 is not enough. Moreover, the proposed SOIGICA algorithm worked well with the additional variability induced by randomized IV. Thus, we recommend performing SOI-GICA by

integrating multiple TC-GICA3 with both different effective SCOs and randomized IVs. The simulated and real fMRI experiments demonstrated that SOIGICA produced more accurate individual- and group-level results than those from traditional TC-GICA3. In particular, SOI-GICA detects group differences with much higher sensitivity and specificity, implying that it is more suitable to detect group differences. As demonstrated by theoretical analyses, the subgroup size (i.e., S) in second-stage PCA reduction is an important parameter in MS-PCA3. As S increases, each subgroup contains more subjects, and the secondstage PCA reduction can use more information to calculate principal components. Thus, MS-PCA3 will produce much better approximated results to MS-PCA1 results. The numerical simulations also demonstrated that TC-GICA3 produced more accurate results with less variability as S increased. However, the SOI-GICA algorithm was still robust against S. For a typical experiment with 30 subjects and S = 4, the total number of different effective SCOs is ∼ 2.7 × 1018. Integrating such a huge number of TC-GICA3 runs is not practical and not necessary. Using numerical simulations, the convergence property of SOI-GICA was evaluated, and an empirical and conservative T value (T = 100) was derived. This T value is quite small when compared with the total number of different effective SCOs, indicating the fast convergence speed of SOI-GICA. Such a relatively fast convergence speed provides a practical guidance for real fMRI applications. Of note, for studies with small group sizes, the problem of SCOinduced variability does not exist because TC-GICA with MS-PCA2 can be adopted, and MS-PCA2 does not introduce SCO-related variability. However, recent TC-GICA studies often involve larger group sizes (e.g., 100 subjects in Stevens et al., 2009). For some multi-center studies, the group size even reaches N 1000 (e.g., 1414 subjects in Biswal et al., 2010). In these cases, MS-PCA3 is a practical choice but is subject to SCO-induced variability. Our proposed SOI-GICA provides a general GICA calculation framework, which provides accurate and reliable results even in very large group studies with limited computational capability. While having many merits, SOI-GICA also has several limitations. First, SOI-GICA integrates results of multiple TC-GICA3 and therefore needs more calculating time (1.5–2 days, mainly depending on group size and component number). Second, in this paper, we only focused on five intensively-studied RSNs that were identified from each TCGICA3 run by using a template-matching approach. If researchers are interested in all estimated components, the problem of component-

1422

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

Fig. 9. Comparison of the estimated spatial maps of the DMN component between SOI-GICA and TC-GICA3 (real resting-state fMRI data set). VIS: the visual component; MOT: sensorimotor component; AUD: auditory component; DMN: default-mode network; and EXE: executive control network. First row: SOI-GICA result (red squares in MDS plots); second and third rows: results from two randomly selected TC-GICA3 runs (circled by green in MDS plots). The yellow arrows indicate the brain areas where TC-GICA3 either failed to detect or mistakenly detected. MDS plots show the relationship between results from SOI-GICA (red squares) and those from multiple TC-GICA3 runs (black dots).

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

matching across different TC-GICA3 runs arises. In our recently released SOI-GICA calculating toolbox “SOI-GICAT” v1.1 beta (http:// www.nitrc.org/projects/cogicat/), we made an initial effort to address each of the above issues. First, the calculating time can be effectively shortened by 40–50% by algorithm optimization (see the manual of SOI-GICAT). Second, an automatic component-matching algorithm was implemented to simultaneously investigate all components without requiring prior templates. Furthermore, the information of subject concatenation order-dependence can be used to rank components with different levels of stability.

Conclusions The results of temporally concatenated group ICA combining more than two stages of PCA reduction depend on subject-concatenation order. Such dependence could reduce the reproducibility of group ICA results between different labs and even impair their reliability. The proposed SOI-GICA approach effectively solved this problem by producing more accurate results than the traditional TC-GICA3, especially in group-difference detection. Moreover, the SOI-GICA algorithm is more robust against various parameters than a single TC-GICA3. Thus, it is highly recommended for large group studies requiring more than two stages of PCA reductions.

2

3 Y11 6 Y12 7 7 6 6 ⋮ 7 7 6 6 Y1K 7 7 6 6 Y21 7 −1 7 = FMSPCA1 6 6 Y22 7; 7 6 6 ⋮ 7 7 6 6Y 7 6 2K 7 4 ⋮ 5 YNK 2

XMSPCA1

(2) where X(1,n) MS − PCA2 denotes the reduced data for subject n, and XMS − PCA2 is the reduced data from the second stage of PCA reduction.

B. Theory on dependency of subject order Here, we provide a theoretical proof of the existence of variability introduced by different subject concatenation orders while performing MS-PCAn (n ≥ 3). Without loss of generality, consider a heuristic example including three random variables (RVs: x1, x2, x3) with zero means and unit variances. We assume correlations 0 b |r(i,j) = corr(xi, xj)| b 1, i ≠ j, i, j = 1, 2, 3. Otherwise, the question becomes trivial. We begin by proving a theorem of PCA reduction on two RVs. Theorem 1. PCA reduction on two RVs Consider two RVs: x1,x2; xg(1,2) denotes the first principal component after PCA reduction of x1 and x2. We have ð1; 2Þ

= rð1;

2Þ

ð1;nÞ

XMSPCA2

x1 + jrð1;

2Þ

jx 2 :

ðB1Þ

Proof 1. Given  Σ = cov ðx1 ; x2 Þ =

1 r ð1;

2Þ

 r ð1; 2Þ ; 1

its eigen-values λ1, λ2 and vectors ω1, ω2 can be calculated based on the normal equation det (Σ − λI) = 0: λ1 = 1 + jrð1; 2Þ ;

2Þ

j;λ2 = 1−jrð1; 2Þ j;

jrð1; 2Þ j

T


; ω2 = jrð1;

2Þ

j;rð1; 2Þ

T

:

Easily, we get the first principal component: xgð1;

2Þ

= ðx1 ; x2 Þω1 = rð1;

2Þ x1

+ jrð1;

2Þ

jx 2 :

Now, we will prove the dependency of grouping concatenation order. Theorem 2. Dependency of grouping ðA1Þ

where F−1 MS − PCA1 is an M × NK reducing matrix, and XMS − PCA1 is the reduced data. MS-PCA2 can be formulized as Eqs. (A2) and (A3): 3 2 Y
−1 6 n1 7 ð1;nÞ 6 Yn2 7; ðn = 1; 2…; NÞ = FMSPCA2 4 ⋮ 5 YnK

ðA3Þ

ð1;NÞ


ω1 = rð1;

Here we provide formulized forms of MS-PCAn (n = 1, 2) procedures. Let Ynk denote brain volume with V voxels observed at the time point k (k = 1, 2, …, K) from subject n (n = 1, 2, …, N). MS-PCA1 can be formulized by:

3

XMSPCA2

xg

Appendix A. Formulization of MS-PCA1 and MS-PCA2

ð1;1Þ

XMSPCA2

7 6 7
−1 6 ð1;2Þ ð2Þ ð2Þ 6 XMSPCA2 7 XMSPCA2 = FMSPCA2 7; 6 7 6 ⋮ 5 4

Acknowledgments This work was supported by the Natural Science Foundation of China (30970773, 30500130, 30621130074 and 30930031), the National Key Basic Research and Development Program (973) Grant No. 2003CB716101 and the Engagement Fund of Outstanding Doctoral Dissertation of Beijing Normal University (No. 08048). We thank Dr. Vince D. Calhoun and three anonymous reviewers for their insightful comments and suggestions. We also thank Mr. Jonathan Adelstein for his assistance in the final language editing.

1423

ðA2Þ

Consider three RVs x1, x2, x3 entered into the second-stage subgrouplevel PCA reduction; there are three different subgrouping methods denoted by ((1, 2), 3) = ((x1, x2), x3), (1, (2, 3)) = (x1, (x2, x3)), and ((1, 3), 2) = ((x1, x3), x2), respectively. We have
2
2
xgðð1;2Þ;3Þ −xgð1;ð2;3ÞÞ + xg ðð1;3Þ;2Þ −xg ð1;ð2;3ÞÞ + xg ðð1;2Þ;3Þ −xg ðð1;

2 3Þ;2Þ

N 0:

ðB2Þ Proof 2. From Theorem 1, we have xgðð1;2Þ;3Þ = rðð1;2Þ;3Þ xgð1;2Þ + jrðð1;2Þ;3Þ jx 3 = rðð1;2Þ;3Þ rð1;2Þ x1 + rðð1;2Þ;3Þ jrð1;2Þ jx 2 + jrðð1;2Þ;3Þ jx 3 :

ðB3Þ ðB4Þ

1424

H. Zhang et al. / NeuroImage 51 (2010) 1414–1424

Similarly, xgð1;ð2;3ÞÞ = rð1;ð2;3ÞÞ x1 + jrð1;ð2;3ÞÞ jrð2;3Þ x2 + jrð1;ð2;

3ÞÞ

jjrð2;3Þ jx 3 ;

xgðð1;3Þ;2Þ = rðð1;3Þ;2Þ rð1;3Þ x1 + jrðð1;3Þ;2Þ jx2 + rðð1;3Þ;2Þ jrð1;3Þ jx3 :

ðB5Þ

ðB6Þ

If the conclusion does not hold, that is, if the PCA result is independent of subject concatenation order (i.e. subgrouping method), then xg((1,2),3) = xg(1,(2,3)) = xg((1,3),2). In other words, rðð1;2Þ;3Þ rð1;2Þ = rð1;ð2;3ÞÞ = rðð1;3Þ;2Þ rð1;3Þ ; rðð1;2Þ;3Þ jrð1;2Þ j = jrð1;ð2;3ÞÞ jrð2;3Þ = jrðð1;3Þ;2Þ j;

jrðð1;2Þ;3Þ j = jrð1;ð2;3ÞÞ jjrð2;3Þ j = rðð1;3Þ;2Þ jrð1;3Þ j: Thus, |r(1,2)| = |r(2,3)| = |r(1,3)| = 1. This result conflicts with our assumption and therefore the conclusion holds, i.e., the existence of a different reduction is proven. Appendix C. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2010.03.039. References Beckmann, C.F., DeLuca, M., Devlin, J.T., Smith, S.M., 2005. Investigations into restingstate connectivity using independent component analysis. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360, 1001–1013. Beckmann, C.F., Mackay, C.E., Filippini, N., Smith, S.M., 2009. Group comparison of resting-state FMRI data using multi-subject ICA and dual regression. NeuroImage 47, S42–S196 Abstract No.441 SU-AM for the Organization for Human Brain Mapping 2009 Annual Meeting. June 18-23, 2009 San Francisco, CA, USA. Biswal, B.B., Mennes, M., Zuo, X.N., Gohel, S., Kelly, C., Smith, S.M., Beckmann, C.F., Adelstein, J.S., Buckner, R.L., Colcombe, S., Dogonowski, A.M., Ernst, M., Fair, D., Hampson, M., Hoptman, M.J., Hyde, J.S., Kiviniemi, V.J., Kotter, R., Li, S.J., Lin, C.P., Lowe, M.J., Mackay, C., Madden, D.J., Madsen, K.H., Margulies, D.S., Mayberg, H.S., McMahon, K., Monk, C.S., Mostofsky, S.H., Nagel, B.J., Pekar, J.J., Peltier, S.J., Petersen, S.E., Riedl, V., Rombouts, S.A., Rypma, B., Schlaggar, B.L., Schmidt, S., Seidler, R.D., Siegle, G.J., Sorg, C., Teng, G.J., Veijola, J., Villringer, A., Walter, M., Wang, L., Weng, X.C., Whitfield-Gabrieli, S., Williamson, P., Windischberger, C., Zang, Y.F., Zhang, H.Y., Castellanos, F.X., Milham, M.P., 2010. Toward discovery science of human brain function. Proc Natl Acad Sci U S A 107, 4734–4739. Calhoun, V.D., Adali, T., Pearlson, G.D., Pekar, J.J., 2001. A method for making group inferences from functional MRI data using independent component analysis. Hum. Brain Mapp. 14, 140–151. Calhoun, V.D., Maciejewski, P.K., Pearlson, G.D., Kiehl, K.A., 2008. Temporal lobe and "default" hemodynamic brain modes discriminate between schizophrenia and bipolar disorder. Hum. Brain Mapp. 29, 1265–1275. Calhoun, V.D., Liu, J., Adali, T., 2009. A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data. NeuroImage 45, S163–S172.

Chen, S., Ross, T.J., Zhan, W., Myers, C.S., Chuang, K.S., Heishman, S.J., Stein, E.A., Yang, Y., 2008. Group independent component analysis reveals consistent resting-state networks across multiple sessions. Brain Res. 1239, 141–151. Damoiseaux, J.S., Rombouts, S.A., Barkhof, F., Scheltens, P., Stam, C.J., Smith, S.M., Beckmann, C.F., 2006. Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U. S. A. 103, 13848–13853. Damoiseaux, J.S., Beckmann, C.F., Arigita, E.J., Barkhof, F., Scheltens, P., Stam, C.J., Smith, S.M., Rombouts, S.A., 2008. Reduced resting-state brain activity in the "default network" in normal aging. Cereb. Cortex 18, 1856–1864. Duann, J.-R., Jung, T.-P., Makeig, S., Sejnowski, T.J., 2003. Consistency of infomax ICA decomposition of functional brain imaging data. Proceedings of the 4th International Symposium on Independent Component Analysis and Blind Signal Separation (ICA 2003). Nara, Japan, pp. 289–294. April. Greicius, M.D., Srivastava, G., Reiss, A.L., Menon, V., 2004. Default-mode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI. Proc. Natl. Acad. Sci. U. S. A. 101, 4637–4642. Grimmett, G.R., Stirzaker, D.R., 1992. Probability and Random Processes, 2nd Edition. Clarendon Press, Oxford. Guo, Y., Pagnoni, G., 2008. A unified framework for group independent component analysis for multi-subject fMRI data. NeuroImage 42, 1078–1093. Himberg, J., Hyvarinen, A., Esposito, F., 2004. Validating the independent components of neuroimaging time series via clustering and visualization. NeuroImage 22, 1214–1222. Hu, D., Yan, L., Liu, Y., Zhou, Z., Friston, K.J., Tan, C., Wu, D., 2005. Unified SPM-ICA for fMRI analysis. NeuroImage 25, 746–755. Hyvarinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural. Netw. 10, 626–634. Jafri, M.J., Pearlson, G.D., Stevens, M., Calhoun, V.D., 2008. A method for functional network connectivity among spatially independent resting-state components in schizophrenia. NeuroImage 39, 1666–1681. Li, Y.O., Adali, T., Calhoun, V.D., 2007. Estimating the number of independent components for functional magnetic resonance imaging data. Hum. Brain Mapp. 28, 1251–1266. Lv, Y.T., Yang, H., Wang, D.Y., Li, S.Y., Han, Y., Zhu, C.Z., He, Y., Tang, H.H., Gong, Q.Y., Zang, Y.F., 2008. Correlations in spontaneous activity and gray matter density between left and right sensorimotor areas of pianists. NeuroReport 19, 631–634. Kiviniemi, V., Kantola, J.H., Jauhiainen, J., Hyvarinen, A., Tervonen, O., 2003. Independent component analysis of nondeterministic fMRI signal sources. NeuroImage 19, 253–260. Kiviniemi, V., Starck, T., Remes, J., Long, X., Nikkinen, J., Haapea, M., Veijola, J., Moilanen, I., Isohanni, M., Zang, Y.F., Tervonen, O., 2009. Functional segmentation of the brain cortex using high model order group PICA. Hum. Brain Mapp. 30, 3865–3886. McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6, 160–188. Mead, A., 1992. Review of the development of multidimensional scaling methods. Statistician 41, 27–39. Smith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E., Johansen-Berg, H., Bannister, P.R., Luca, M.D., Drobnjak, I., Flitney, D.E., Niazy, R., Saunders, J., Vickers, J., Zhang, Y., Stefano, N.D., Brady, J.M., Matthews, P.M., 2004. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage 23, 208–219. Smith, S.M., Fox, P.T., Miller, K.L., Glahn, D.C., Fox, P.M., Mackay, C.E., Filippini, N., Watkins, K.E., Toro, R., Laird, A.R., Beckmann, C.F., 2009. Correspondence of the brain's functional architecture during activation and rest. Proc. Natl. Acad. Sci. U. S. A. 106, 13040–13045. Stevens, M.C., Pearlson, G.D., Calhoun, V.D., 2009. Changes in the interaction of resting-state neural networks from adolescence to adulthood. Hum. Brain Mapp. 30, 2356–2366. Yang, Z., Laconte, S., Weng, X., Hu, X., 2008. Ranking and averaging independent component analysis by reproducibility (RAICAR). Hum. Brain Mapp. 29, 711–725. Ylipaavalniemi, J., Vigario, R., 2008. Analyzing consistency of independent components: an fMRI illustration. NeuroImage 39, 169–180. Zuo, X.N., Kelly, C., Adelstein, J.S., Klein, D.F., Castellanos, F.X., Milham, M.P., 2010. Reliable intrinsic connectivity networks: test–retest evaluation using ICA and dual regression approach. NeuroImage 49, 2163–2177.