A Simple and Objective Method for Reproducible Resting State Network (RSN) Detection in fMRI

Spatial Independent Component Analysis (ICA) decomposes the time by space functional MRI (fMRI) matrix into a set of 1-D basis time courses and their associated 3-D spatial maps that are optimized for mutual independence. When applied to resting state fMRI (rsfMRI), ICA produces several spatial independent components (ICs) that seem to have biological relevance - the so-called resting state networks (RSNs). The ICA problem is well posed when the true data generating process follows a linear mixture of ICs model in terms of the identifiability of the mixing matrix. However, the contrast function used for promoting mutual independence in ICA is dependent on the finite amount of observed data and is potentially non-convex with multiple local minima. Hence, each run of ICA could produce potentially different IC estimates even for the same data. One technique to deal with this run-to-run variability of ICA was proposed by [1] in their algorithm RAICAR which allows for the selection of only those ICs that have a high run-to-run reproducibility. We propose an enhancement to the original RAICAR algorithm that enables us to assign reproducibility -values to each IC and allows for an objective assessment of both within subject and across subjects reproducibility. We call the resulting algorithm RAICAR-N (N stands for null hypothesis test), and we have applied it to publicly available human rsfMRI data (http://www.nitrc.org). Our reproducibility analyses indicated that many of the published RSNs in rsfMRI literature are highly reproducible. However, we found several other RSNs that are highly reproducible but not frequently listed in the literature.


Introduction
Independent component analysis (ICA) [2][3][4][5] models the observed data as a linear combination of a set of statistically independent and unobservable sources [6]. first proposed the application of ICA to the analysis of functional magnetic resonance imaging (fMRI) data. Subsequently, ICA has been applied to fMRI both as an exploratory tool for the purpose of identifying task related components [6] as well as a signal clean up tool for the purpose of removing artifacts from the fMRI data [7]. Recently, it has been shown that ICA applied to resting state fMRI (rsfMRI) in healthy subjects reveals a set of biologically meaningful spatial maps of independent components (ICs) that are consistent across subjectsthe so called resting state networks (RSNs) [8]. Hence, there is a considerable interest in applying ICA to rsfMRI data in order to define the set of RSNs that characterize a particular group of human subjects, a disease, or a pharmacological effect.
Several variants of the linear ICA model have been applied to fMRI data including square ICA (with equal number of sources and sensors) [9], non-square ICA (with more sensors than sources) [6], and non-square ICA with additive Gaussian noise (noisy ICA) [10]. All of these models are well known in the ICA literature [2,3,5,11]. Since the other ICA models are specializations of the noisy ICA model, we will assume a noisy ICA model henceforth.
Remarkably, the ICA estimation problem is well posed in terms of the identifiability of the mixing matrix given several non-Gaussian and at most 1 Gaussian source in the overall linear mixture [3,[12][13][14]. In the presence of more than 1 Gaussian source, such as in noisy ICA, the mixing matrix corresponding to the non-Gaussian part of the linear mixture is identifiable (upto permutation and scaling). In addition, the source distributions are uniquely identifiable (upto permutation and scaling) given a noisy ICA model with a particular Gaussian co-variance structure, for example, the isotropic diagonal co-variance. For details, see section 2.1.2.
While these uniqueness results are reassuring, a number of practical difficulties prevent the reliable estimation of ICs on real data. These difficulties include (1) true data not describable by an ICA model, (2) ICA contrast function approximations, (3) multiple local minima in the ICA contrast function, (4) confounding Gaussian noise and (5) model order overestimation. See section 2.1.3 for more details. A consequence of these difficulties is that multiple ICA runs on the same data or different subsets of the data produce different estimates of the IC realizations.
One technique to account for this run-to-run variability in ICA was proposed by [15] in their algorithm ICASSO. Using repeated runs of ICA with bootstrapped data using various initial conditions, ICASSO clusters ICs across ICA runs using agglomerative hierarchical clustering and also helps in visualizing the estimated ICs. The logic is that reliable ICs will show up in almost all ICA runs and thus will form a tight cluster well separated from the rest. [16] proposed a technique similar to ICASSO called selforganizing group ICA (sogICA) which allows for clustering of ICs via hierarchical clustering in across subject ICA runs. When applied to multiple ICA runs across subjects, ICASSO does not restrict the IC clusters to contain only 1 IC from each subject per ICA run. In contrast, sogICA allows the user to select the minimum number of subjects for a ''group representative'' IC cluster containing distinct subjects. By labelling each ICA run as a different ''subject'' sogICA can also be applied to analyze multiple ICA runs across subjects.
Similar in spirit to ICASSO and sogICA [1], proposed an intuitive approach called RAICAR (Ranking and Averaging Independent Component Analysis by Reproducibility) for reproducibility analysis of estimated ICs. The basic idea in RAICAR is to select only those ICs as ''interesting'' or ''stable'' which show a high run-to-run ''reproducibility''. RAICAR uses simple and automated spatial cross-correlation matrix based IC alignment, which has been shown to be more accurate compared to ICASSO [1]. RAICAR is applicable to both within subject as well as across subjects reproducibility analysis.
A few limitations of ICASSO, sogICA and RAICAR are worth noting: N ICASSO requires the user to select the number of IC clusters and is inapplicable without modification for across subjects analysis of ICA runs since the IC clusters are not restricted to contain only 1 IC per ICA run.
N sogICA requires the user to select the minimum number of subjects for a ''group representative'' cluster and also a cutoff on within cluster distances.
N RAICAR uses an arbitrary threshold on the reproducibility indices selected ''by eye'' or set at an arbitrary value, such as 50% of the maximum reproducibility value.
We propose a simple extension to RAICAR that avoids subjective user decisions and allows for an automatic reproducibility cutoff. The reproducibility indices calculated in RAICAR differ in magnitude significantly depending on whether the input to RAICAR: N (a) is generated using multiple ICA runs on the same data N (b) comes from multiple ICA runs on varying data sets (e.g. between and across subject runs) See Figure 1 for an illustration of this effect. Obviously, the reproducibility indices are much lower in case (b) since we account for both within subject and between subjects variability in estimating ICs. Case (b) is also of great interest from a practical point of view since we are often interested in making statements about a group of subjects. Hence, it is clear that a cutoff on RAICAR reproducibility values for the purposes of selecting the ''highly reproducible'' components should be data dependent. In this work, 1. We propose a modification of the original RAICAR algorithm by introducing an explicit ''null'' model of no reproducibility. 2. We use this ''null'' model to automatically generate p-values for each IC via simulation. This allows for an objective cutoff specification for extracting reproducible ICs (e.g. reproducible at pv0:05) within and across subjects. We call the resulting algorithm RAICAR-N (N stands for ''null'' hypothesis test). 3. We validate RAICAR-N by applying it to publicly available human rsfMRI data.

Notation
N The set of real numbers will be denoted by R. Scalars variables and functions will be denoted in a non-bold font (e.g., s 2 ,L,p or Y,f ). Vectors will be denoted in a bold font (except Greek letters) using lower case letters (e.g., y,m,g). Matrices will be denoted in bold font using upper case letters (e.g., A,S,W). The transpose of a matrix A will be denoted by A T and its inverse will be denoted by A {1 . I p will denote the p|p identity matrix and 0 will denote a vector or matrix of all zeros whose size should be clear from context. N L is the number of ways of choosing L objects from N objects when order does not matter.
N The jth component of vector t i will be denoted by t ij whereas the jth component of vector t will be denoted by t j . The element (i,j) of matrix G will be denoted by G(i,j) or G ij .
Estimates of variables will be denoted by putting a hat on top of the variable symbol. For example, an estimate of s will be denoted byŝ s.
N If x is a random vector with a multivariate Normal distribution with mean m and covariance S then we will denote this distribution by N xjm,S ð Þ. The joint density of vector s will be denoted by p s (s) whereas the marginal density of s i will be denoted as p si (s i ). E f (s,h) ½ denotes the expectation of f (s,g) with respect to both random variables s and g.

Methods
The organization of this article revolves around the following sequence of questions, which ultimately lead to the development of RAICAR-N: 1. Why is a reproducibility assessment necessary in ICA analysis?
In order to answer this question, we cover the fundamentals of ICA including identifiability issues in sections 2.

ICA background
In this section, we provide a brief introduction to ICA along with a discussion of associated issues related to model order selection, identifiability and run-to-run variability. The noisy ICA model assumes that observed data y is generated as a linear combination of unobservable independent sources confounded with Gaussian noise: In this model, If the marginal density of the ith source s i is p si (s i ) then the joint source density p s (s) factorizes as P q i~1 p si (s i ) because of the independence assumption but is otherwise assumed to be unknown. Also, since the elements of s are independent their co-variance matrix D is diagonal. The set of variables F~m,A,D,S f grepresents the unknown parameters in the noisy ICA model. Before discussing the identifiability of model 2.1, we briefly discuss the choice of model order or the assumed number of ICs q.
2.1.1 Estimating the model order q. Rigorous estimation of the model order q in noisy ICA is difficult as the IC densities p si (s i ) are unknown. This means that p yjq,F ð Þ, the marginal density of the observed data given the model order and the ICA parameters cannot be derived in closed form (by integrating out the ICs) without making additional assumptions on the form of IC densities. Consequently, standard model selection criteria such as Bayes information criterion (BIC) [17] cannot be easily applied to the noisy ICA model to estimate q. One solution is to use a factorial mixture of Gaussians (MOG) joint source density model as in [5], and use the analytical expression for p yjq,F ð Þ in conjunction with BIC. This solution is quite general in terms of allowing for an arbitrary Gaussian noise co-variance S, but maximizing p yjq,F ð Þwith respect to F becomes computationally intractable using an expectation maximization (EM) algorithm for qw13 ICs [5]. Another rigorous non-parametric approach for estimating q that is applicable to the noisy ICA model with isotropic diagonal Gaussian noise co-variance i.e., with S~s 2 I p is the random matrix theory based sequential hypothesis testing approach of [18]. To the best of our knowledge, these are the only 2 rigorous approaches for estimating q in the noisy ICA model.
Approximate approaches for estimating q commonly used in fMRI literature (e.g., [10]) consist of first relaxing the isotropic diagonal noisy ICA model (with S~s 2 I p ) into a probabilistic PCA (PPCA) model of [19] where the source densities are assumed to be Gaussian i.e., where p s (s)~N sj0,I q À Á . When using the PPCA model, it becomes possible to integrate out the Gaussian sources to get an expression for p yjq,F ð Þthat can be analytically maximized [19]. Subsequently, methods such as BIC can be applied to estimate q. Alternative approaches for estimating q in the PPCA  [20], or in datarich situations such as fMRI, even the standard technique of crossvalidation [21].
From a biological point of view, it has been argued [22] that the number of extracted ICs simply reflect the various equally valid views of the human functional neurobiology -smaller number of ICs represent a coarse view while a larger number of ICs represent a more fine grained view. However, it is worth noting that from a statistical point of view, over-specification of q will lead to overfitting of the ICA model, which might render the estimated ICs less generalizable across subjects. On the other hand, underspecification of q will result in incomplete IC separation. Both of these scenarios are undesirable.
2 Thus the mean vector m is exactly identifiable. Identifiability of A: A fundamental decomposition result states that the noisy ICA problem is well-posed in terms of the identifiability of the mixing matrix A upto permutation and scaling provided that the components of s are independent and non-Gaussian [3,[12][13][14]. If L is a diagonal scaling matrix and P is a permutation matrix then the identifiability result can be stated as: where 2.3 is another decomposition of y with s 1 containing independent and non-Gaussian components. In other words, the mixing matrix A is identifiable upto permutation and scaling. Identifiability of D and S: Equating the second moments of the right hand side of 2.3 and 2.1 and noting the equality of means 2.5 and the independence of s,h and s 1 ,h 1 we get: Let W be a q|p matrix andQ Q be a p|(p{q) orthogonal matrix such that: From 2.8 and 2.7 we get: Case 1: S~s 2 I p and S 1~s 2 1 I p . The second equation in 2.9 along with the orthogonality ofQ Q gives s 2~s2 1 and thus S~S 1 . If we fix the scaling of A 1 by selecting L 2~I q then from the first equation in 2.9 we get:

ð2:10Þ
In other words, the noise co-variance S~s 2 I p is uniquely determined and for a fixed scaling L 2~I q , the source variances D are also uniquely determined upto permutation. Case 2: S and S 1 arbitrary positive definite matrices. Suppose X is a square matrix and let diag(X) be the diagonal matrix obtained by setting the non-diagonal elements of X to 0 and similarly let offdiag(X) be the matrix obtained by setting the diagonal elements of X to 0. The noise-covariance is partially identifiable by the following conditions: For a fixed scaling L 2~I q , the sources variances D,D 1 are constrained by: In general, the source variances D cannot be uniquely determined as noted in [14]. Identifiability of the distribution of s: Is the distribution of the non-Gaussian components of s identifiable? From 2.1 and 2.3: Noting the independence of s,g and s 1 ,g 1 : Now g and g 1 are multivariate Gaussian random vectors both with mean 0 and co-variance matrix S and S 1 respectively. Hence, their characteristic functions are given by [23,24]: Claim 2.1 A sufficient condition for identifiability upto permutation and scaling of the non-Gaussian distributions in s given two different parameterizations in 2.1 and 2.3 is: ProofN From 2.20 and 2.11, we get: Thus from 2.19, are not equal to 0 for any finite t, therefore, from 2.22 and 2.18 we get: Note that L is a diagonal scaling matrix with entries l 1 ,l 2 , . . . ,l q on the diagonal and P is a permutation matrix. Thus, . .
where i 1 ,i 2 , . . . ,i q is some permutation of integers 1,2, . . . ,q. Suppose Y s(j) is the characteristic function of the jth component of s and Y s 1 (j) is the characteristic function of the jth component of s 1 . Since the components of s and s 1 are independent by assumption, the joint characteristic functions Y(s) and Y(s 1 ) factorize: From 2.25 and 2.23 All characteristic functions satisfy [23,24]: Since i 1 ,i 2 , . . . ,i q is simply a permutation of integers 1,2, . . . ,q, there exists a j such that i j~1 . Then set t 2~0 ,t 3~0 , . . . ,t q~0 in 2.26. Then 2.27 and 2.26 imply: Select the scaling matrix as L 2~I q and thus L is a diagonal matrix with elements +1 on the diagonal. Thus l 1~+ 1 and 2.28 can be re-written as: Therefore, Hence the characteristic function of the 1st component of s is identical to the characteristic function of the (possibly sign-flipped) jth component of s 1 . Since characteristic functions uniquely characterize a probability distribution [23], the distribution of s(1) and +s 1 (j) is identical. Next, by setting t 1~0 ,t 3~0 , . . . ,t q~0 , we can find a distribution from s 1 that matches the 2nd component  While the source distributions might not be uniquely identifiable for arbitrary co-variance matrices S, they are indeed uniquely identifiable upto permutation and scaling for the noisy ICA model with isotropic Gaussian noise co-variance. For more general conditions that guarantee uniqueness of source distributions, please see [25,26].
Corollary 2.2 If S~s 2 I p and S 1~s 2 1 I p , then the source distributions are uniquely identifiable upto sign flips for L 2~I q . Proof Suppose S~s 2 I p and S 1~s 2 1 I p . Then from 2.9 S~S 1 and thus diag(WSW T )~diag(WS 1 W T ). The corollary then follows from Claim 2.1.
Corollary 2.3 If D~D 1~Iq , then the source distributions are uniquely identifiable up to sign flips for L 2~I q . Proof If D~D 1~Iq , then noting that PP T~I q , we get D~PD 1 P T . Hence from 2.12, we get diag(WSW T )~diag (WS 1 W T ). The corollary then follows from Claim 2.1.

Why is there a run-to-run variability in estimated
ICs?. From the discussion in section 2.1.2, it is clear that for a noisy ICA model with isotropic diagonal additive Gaussian noise co-variance: 1. The noisy ICA parameters F~m,A,D,S f g are uniquely identifiable up to permutation and scaling. 2. The source distributions in s are uniquely identifiable upto permutation and scaling.
While the above theoretical properties of ICA are reassuring, there are a number of practical difficulties that prevent the reliable estimation of ICs on real data: 1. Validity of the ICA model: The assumption that the observed real data is generated by an ICA model is only that -an ''assumption''. If this assumption is not valid, then the uniqueness results do not hold anymore. 2. Mutual information approximations: From an information theoretic point of view, the ICA problem is solved by minimizing a contrast function which is an approximation to the mutual information [27] between the ICs that depends on the finite amount of observed data. Such an approximation is necessary, since we do not have access to the marginal source densities p si . Different approximations to mutual information will lead to different objective functions and hence different solutions. This is one of the reasons why different ICA algorithms often produce different IC estimates even for the same data. 3. Non-convexity of ICA objective functions: The ICA contrast function is potentially non-convex and hence has multiple local minima. Since global minimization is a challenging problem by itself, most ICA algorithms will only converge to local minima of the ICA contrast function. The run-to-run variability of IC estimates will also depend on the number of local minima in a particular ICA contrast function. 4. IC estimate corruption by Gaussian noise: For noisy ICA, the IC realizations cannot be recovered exactly even if the true mixing matrix A and mean vector m are known in 2.1. Commonly used estimators for recovering realization of ICs include the least squares [10] as well as the minimum mean square error (MMSE) [14]. Consider the least squares estimateŝ s of a realization of s based on y: This means that even for known parameters, IC realization estimatesŝ s will be corrupted by correlated Gaussian noise.
Hence using different subsets of the data under the true model will also lead to variability in estimated ICs. 5. Over-fitting of the ICA model: Over specification of the model order leads to the problem of over-fitting in ICA. As we describe below, this can lead to (1) the phenomenon of IC ''splitting'' and (2) an increase in the variance of the IC estimates.

IC ''splitting''
Suppose that the true model order or the number of non-Gaussian sources in an ICA decomposition of y such as 2.1 is q. Then a fundamental result in [12,Theorem 1] states that for any other ICA decomposition of y, the number of non-Gaussian sources remains the same while the number of Gaussian sources can change. In other words, y cannot have two different ICA decompositions containing different number of non-Gaussian sources.
In view of this fact, how can a model order q ICA decomposition containing q non-Gaussian sources be ''split'' into a (qz1) ICA decomposition containing (qz1) non-Gaussian sources when performing ICA estimation using an assumed model order of (qz1)? As we describe below, the order (qz1) ICA decomposition is only an approximation to the order q ICA decomposition.
Let a i be the ith column of A in 2.1. In the presence of noise, it might be possible to approximate: Here:   i . By replacing a i s i in 2.1 using 2.32, we arrive at an approximate model order (qz1) decomposition of y. In this decomposition, the component s i from a model order q decomposition appears to be ''split'' into two sub-components: s 1 i and s 2 i .

Inflated variance of IC estimates
Overestimation of model order will lead to over-fitting of the mixing matrix A. In other words, A could have several columns that are highly correlated with each other. This could happen as a result of IC ''splitting'' as discussed above. Now, for a given realization s, the variance ofŝ s is given by Var(ŝ s)~s 2 (A T A) {1 (for isotropic Gaussian co-variance). An increase in number of columns of A and the fact that many of them are highly correlated implies that the variability of IC estimates Var(ŝ s) is inflated.
In other words, running ICA multiple times on the same data or variations thereof with random initialization could produce different ICs.

ICA algorithms, single subject ICA and group ICA
In this section, we give a brief summary of how the ICA parameters are estimated in practice and also summarize the two most common modes of ICA application to fMRI data -single subject ICA (section 2.2.1) and temporal concatenation based group ICA (section 2.2.2).
Given several independent observations y as per the noisy ICA model 2.1, most ICA algorithms estimate the ICA parameters F~m,A,D,S f g and the realizations of s in 2 steps. We only consider the case with S~s 2 I p , since as shown in section 2.1.2, the mixing matrix A and source distributions of s are identifiable upto permutation and scaling for this case.
1. First, the diagonal source co-variance is arbitrarily set as D~I q . The mean vector m is estimated as E y ð Þ. Then, using PCA or PPCA [19], the mixing matrix A is estimated, upto an orthogonal rotation matrix O, to be in a signal subspace which is spanned by the principal eigenvectors corresponding to the largest eigenvalues of the data co-variance matrix E (y{m)(y{m) T Â Ã . The noise variance s 2 is estimated in this step as well. 2. Next, an estimatorŝ s for the source realizations is defined using techniques such as least squares or MMSE. The only unknown involved in these estimates is the orthogonal rotation matrix O. 3. Finally, the non-Gaussianity of the empirical density of components ofŝ s is optimized with respect to O using algorithms such as fixed point ICA [27,28].
For more details on noisy ICA estimation, please see [10] and for more details on ICA algorithms, please see [29].
2.2.1 Single subject ICA. How is ICA applied to single subject fMRI data? Suppose we are given a single subject fMRI scan which we rearrange as a p|n 2D matrix Y in which column i is the p|1 observed time-course y i in the brain at voxel i. Observed time-courses y 1 ,y 2 , . . . ,y n are considered to be n independent realizations of y as per the linear ICA model 2.1. SupposeŜ S~½ŝ s 1 ,ŝ s 2 , . . . ,ŝ s n is the q|n matrix containing the estimated source realizations at the n voxels. The j th row ofŜ S is the jth IC. In other words, we decompose the time by space fMRI 2D matrix into a set of basis time-courses and a set of q 3D IC maps using ICA.
2.2.2 Group ICA. How is ICA applied to data from a group of subjects in fMRI? Suppose we collect fMRI images from m subjects. First, we register all subjects to a common space using a registration algorithm (e.g., affine registration). Next, we rearrange each of the fMRI scans into m 2D matrices Y 1 . . . Y m , each of size p|n. Column j in Y i is the demeaned time-course observed at voxel location j for subject i. The matrices Y 1 . . . Y m are temporally concatenated to get a pm|n matrix Z as follows: Column i of Z is the pm|1 vector z i which is assumed to follow a linear ICA model 2.1. z 1 ,z 2 , . . . ,z n are considered to be independent realizations of the model 2.1. SupposeŜ S G1ŝ s 1 ,ŝ s 2 , . . . ,ŝ s n is a q|n matrix containing the estimated source realizations at the n voxels. The jth row ofŜ S G is the j th group IC. In group ICA, the joined time-series across subjects is modeled using noisy linear ICA. In practice, Y i is the PCA reduced data set for subject i. The PCA reduction is either done separately for each subject using subject specific data co-variance [9] or an average data co-variance across subjects [8]. The average co-variance approach requires each subject to have the same number of time points in fMRI scans.

The original RAICAR algorithm
In this section, we give a brief introduction to the RAICAR algorithm of [1]. Suppose we are given a data set which we decompose into n C ICs using ICA (e.g., single subject or group ICA). Our goal is to assess which ICs consistently show up in multiple ICA runs i.e., the reproducibility of each of these n C ICs. To that extent, we run the ICA algorithm K times. Suppose x (m) j is the n|1 vector (e.g. spatial ICA map re-arranged into a vector) of the j th IC from m th ICA run. Suppose G lm is a n C |n C absolute spatial cross-correlation coefficient matrix between the ICs from runs l and m: where j:j denotes absolute value. G lm (i,j) is the absolute spatial cross-correlation coefficient between IC i from run l and IC j from run m. The matrices G lm are then arranged as elements of a K|K block-matrix G such that the l th row and m th column of G is G(l,m)~G lm (Figure 2). This block matrix G is the starting point for a RAICAR across-run component matching process.
Since ICs within a particular run cannot be matched to each other, the n C |n C matrices G(l,l),l~1 . . . K along the blockdiagonal of G are set to 0 as shown in Figure 2 with a gray color. The following steps are involved in a RAICAR analysis: Similarly, suppose element (i,b s ) is the maximal element in the i th row of G ls . Then component b s is the best matching component from run s with component i from run l. As noted in [1], in most cases a s~bs . However, it is possible that a s =b s . Hence the component number e s matching MC 1 from run s is defined as follows: We would also like to remove component e s of run s from further consideration during the matching process. To that extent, we zero out the e s th row from G sr ,r~1 . . . K and the e s th column from G rs ,r~1 . . . K.

Once a matching component e s has been found for all runs s=l,m,
we also zero out the ith row from G lr ,r~1 . . . K and the ith column from G rl ,r~1 . . . K. Similarly, we zero out the j th column from G rm ,r~1 . . . K and the j th row from G mr ,r~1 . . .
The normalized reproducibility of MC s is then defined as: The double sum in 2.37 is simply the sum of the upper triangular part of H MCs excluding the diagonal. The normalizing factor (K {1)K 2 is simply the maximum possible value of this sum.
Note that our definition of normalized reproducibility is slightly different from that in [1]. Whereas [1] averages the thresholded absolute correlation coefficients, we simply average the unthresholded absolute correlation coefficients to compute reproducibility thereby avoiding the selection of a threshold on the absolute correlation coefficients.

The RAICAR-N enhancement
In this section, we describe how to compute reproducibility pvalues for each matched component in RAICAR. Note that the RAICAR ''component matching'' process can be used to assess the reproducibility of any spatial component maps -not necessarily ICA maps. For instance, RAICAR can be used to assess the reproducibility of a set of PCA maps across subjects.  [1]. The ICA algorithm is run K times with each run producing n C ICs. G is a K|K block matrix with elements G(l,m)~G lm where G lm is the n C |n C absolute spatial cross-correlation matrix between ICs from runs l and m. The numbered green circles indicate the sequence of steps in applying RAICAR to a given data set. Our definition of normalized reproducibility in box 7 averages un-thresholded correlation coefficients thereby avoiding the selection of a correlation coefficient threshold prior to averaging. doi:10.1371/journal.pone.0027594.g002 In order to generate reproducibility p-values for the matched component maps: 1. We need to determine the distribution of normalized reproducibility that we get from the RAICAR ''component matching'' process when the input to RAICAR represents a set of ''non-reproducible component maps'' across the K runs. 2. In addition, we would also like to preserve the overall structure seen in the observed sets of spatial component maps across the K runs when generating sets of ''non-reproducible component maps'' across the K runs.
Hence for IC reproducibility assessment, we propose to use the original set of ICs across the K runs to generate the ''nonreproducible component maps'' across the K runs.
Suppose K ICA runs are submitted to RAICAR which gives us a n C |1 vector of observed normalized reproducibility values Reproducibility(MC i ),i~1 . . . n C -one for each IC. We propose to attach p-values for measuring the reproducibility of each IC in a data-driven fashion as follows: 1. First, we label the Kn C ICs across the K runs using unique integers. In run 1, the ICs are labelled using integers 1, . . . ,n C . In run 2, the ICs are labelled using integers (n C z1), . . . ,2n C and so on. In run K, the ICs are labelled using integers (K{1)n C z1, . . . ,Kn C .

Our ''null'' hypothesis is:
H 0 : None of the ICs are reproducible ð2:38Þ To do this, we randomly permute the integers 1,2, . . . ,Kn C to get the permuted integers p(1), . . . ,p(Kn C ). Obviously p(i)=p(j) if i=j. vector of normalized reproducibility Reproducibility Null under H 0 . 6. Finally, we assign a p-value for reproducibility to each matched IC across the K runs. The observed reproducibility for i th matched IC is Reproducibility(MC i ) and its p-value is: ð2:39Þ 7. Only those components with Reproducibility pval (MC i )vp crit are considered to be significantly reproducible. We can use a fixed and objective value for p crit such as 0:05. Note that this fixed cutoff is independent of the amount of variability in the input to RAICAR-N. Please see Figure 3 for a pictorial depiction of this process. Figure 4 shows a schematic of the RAICAR-N analysis process.

How many subjects should be used per group ICA run in RAICAR-N?
The input to RAICAR-N can either be single subject ICA runs or group ICA runs across a set of subjects. Note that the individual subject ICA runs are spatially unconstrained whereas a group ICA spatially constrains the group ICs across a set of subjects. Hence the number of ICs that can be declared as significantly reproducible at the group level are usually more than those that can be declared significantly reproducible at the single subject level. Hence the following question is relevant: Suppose we have a group of N subjects. We randomly select L subjects and form a single group of subjects. We repeat this process K times to get K groups of L subjects each of which is subjected to a group ICA analysis. Given the number of subjects N, how should we choose L and K?
First, we discuss the choice of L. If L~N then each of the K groups will contain the same N subjects and hence there will be no diversity in the K groups. We would like to control the amount of diversity in the K groups of L subjects. Consider any 2 subjects X and Y . The probability P XY (L) that both X and Y appear in a set of L randomly chosen subjects from N subjects is given by: The expected number of times that X and Y appear together in sets of L subjects out of K independently drawn sets is: Ideally, we would like E XY (L) to be only a small fraction of K. Hence we impose the restriction: where a max is a user defined constant such as a max~0 :05. This implies that the chosen value of L must satisfy: In practice, we choose the largest value of L that satisfies this inequality. As shown in Figure 5, if N~23 and a max~0 :05 then the largest value of L that satisfies 2.43 is L~5. The number of group ICA runs K should be as large as possible. From our experiments on real fMRI data we can roughly say that values of Kw50 give equivalent results.
2.6 How to display the estimated non-Gaussian spatial structure in ICA maps?
The ICs have been optimized for non-Gaussianity. However, there can be many types of non-Gaussian distributions. It has been empirically found that the non-Gaussian distributions of ICs found in fMRI data have the following structure:

A central Gaussian looking part and 2. A tail that extends out on either end of the Gaussian
It has been suggested in [10] that a Gaussian/Gamma mixture model can be fitted to this distribution and the Gamma components can be thought of as representatives of the non-Gaussian structure. We follow a similar approach: 1. The output of a RAICAR-N analysis is a set of spatial ICA maps (either z-transformed maps or raw maps) concatenated into a 4-D volume. 2. We do a voxelwise transformation to Normality using the voxelwise empirical cumulative distribution function as described in [30].
3. Next, we submit the resulting 4-D volume to a voxelwise group analysis using ordinary least squares. The design matrix for group analysis depends on the question being considered. In our case, the design matrix was simply a single group average design. 4. The resulting t-statistic maps are subjected to Student t, Gamma pos and Gamma neg mixture modeling. The logic is that if the original ICA maps are pure Gaussian (i.e., have no interesting non-Gaussian structure) then the result of a group average analysis will be a pure Student t map which will be captured by a single Student t (i.e., the Gamma pos and Gamma neg will be driven to 0 class fractions). Hence the ''null'' hypothesis will be correctly accounted for. . Flowchart for a group ICA based RAICAR-N analysis. The N single subject data sets are first pre-processed and subsequently bootstrapped to create K groups, each group containing L distinct subjects. Each group of L subjects is submitted to a temporal concatenation group ICA analysis. The resulting IC maps (either raw ICs or ICs scaled by noise standard deviation) are subjected to a RAICAR analysis. The crossrealization cross correlation matrix (CRCM) is randomly permuted multiple times: G?G(g,g) where g is a random permutation of integers from 1, . . . ,Kn C . The permuted CRCMs are subjected to a RAICAR analysis to generate a realization of reproducibility values under the ''null'' hypothesis. The computed ''null'' distribution of reproducibility values is used to assign p values to the observed reproducibility of the original RAICAR run. Finally, reproducible ICs are averaged using a random effects analysis and the resulting t-statistic images are subjected to Gamma neg , Student t and Gamma pos mixture modeling. doi:10.1371/journal.pone.0027594.g004 5. If the Gamma distributions have w0:5 posterior probability at some voxels then those voxels are displayed in color to indicate the presence of significant non-Gaussian structure over and above the background Student t distribution.
Examples of Student t, Gamma pos and Gamma neg mixture model fits are shown in Figure 6.

Preprocessing
Data was analyzed using tools from the FMRIB software library (FSL: http://www.fmrib.ox.ac.uk/fsl/). Preprocessing steps included motion correction, brain extraction, spatial smoothing with an isotropic Gaussian kernel of 5 mm FWHM and 100 s high-pass temporal filtering. Spatial ICA was performed using a noisy ICA model as implemented in FSL MELODIC [10] in either single subject or multi-subject temporal concatenation mode also called group ICA. Please see section 2.2 for a brief summary of single subject ICA and group ICA. In each case, we fixed the model order of ICA at q~40 to be consistent with the model order range typically extracted in rsfMRI and fMRI [16,31]. For temporal concatenation based group ICA, single subject data was first affinely registered to the MNI 152 brain and subsequently resampled to 46464 resolution (MNI 46464) to decrease computational load. Please see Figure 4 for a schematic of the RAICAR-N analysis process. In this work, we report across subject RAICAR-N analyses, but as shown in Figure 7, within subject ICA runs can also be entered into RAICAR-N.

RAICAR-N analysis with 1 ICA run per subject
Spatial ICA was run once for each of the N~23 subjects in their native space. The resulting set of ICA components across subjects were transformed to MNI 46464 space and were submitted to a RAICAR-N analysis. In all RAICAR-N analyses reported in this article, we used the z-transformed IC maps -which are basically the raw IC maps divided by a voxelwise estimate of noise standard deviation (named as melodic_IC.nii.gz in ME-LODIC). It is also possible to use the raw IC maps as inputs to RAICAR-N. ICA components were sorted according to their reproducibility and p-values were computed for each ICA component. Please see Figure 8.
We compared the reproducible RSNs from the single subject RAICAR-N analysis to the group RSN maps reported in literature [8]. Please see Figure 9.
To summarize, when single subject ICA runs are combined across subjects:

RAICAR-N on random sets of 5 subjects -50 group ICA runs
To promote diversity across the group ICA runs, as discussed in section 2.5, L~5 subjects were drawn at random from the group of N~23 subjects and submitted to a temporal concatenation based group ICA. This process was repeated K~50 times and the resulting set of 50 group ICA maps were submitted to a RAICAR-N analysis. ICA components were sorted according to their reproducibility and p-values were computed for each ICA component. Please see Figure 10.
We compared the reproducible RSNs from the single subject RAICAR-N analysis to the RSN maps reported in literature [8]. Please see Figure 11.
In summary, when 50 random 5 subject group ICA runs (from a population of 23 subjects) are combined using RAICAR-N:

RAICAR-N on random sets of 5 subjects -100 group ICA runs
To promote diversity across the group ICA runs, as discussed in section 2.5, L~5 subjects were drawn at random from the group of N~23 subjects and submitted to a temporal concatenation based group ICA. This process was repeated K~100 times and the resulting set of 100 group ICA maps were submitted to a RAICAR-N analysis. ICA components were sorted according to  their reproducibility and p-values were computed for each ICA component. Please see Figure 12.
We compared the reproducible RSNs from the single subject RAICAR-N analysis to the RSN maps reported in literature [8]. Please see Figure 13.
In summary, when 100 random 5 subject group ICA runs (from a population of 23 subjects) are combined using RAICAR-N: N There is 1 other ''non-standard'' RSN that achieves a p-value of 0.05824.

Discussion
As discussed in section 2.1.2, in the noisy linear ICA model with isotropic diagonal Gaussian noise co-variance, for a given true model order, the mixing matrix and the source distributions are identifiable upto permutation and scaling. However, as pointed out in section 2.1.3, various factors prevent the convergence of ICA algorithms to unique IC estimates. These factors include ICA model not being the true data generating model, approximations to mutual information used in ICA algorithms, multiple local minima in ICA contrast functions, confounding Gaussian noise as well as variability due to model order over-estimation. A practical implication of these factors is that ICA algorithms converge to different IC estimates depending on how they are initialized and on the specific data used as input to ICA. Hence, there is a need for a rigorous assessment of reproducibility or generalizability of IC estimates. A set of reproducible ICs can then be used as ICA based characteristics of a particular group of subjects.
We proposed an extension to the original RAICAR algorithm for reproducibility assessment of ICs within or across subjects (Figure 7). The modified algorithm called RAICAR-N builds up a ''null'' distribution of normalized reproducibility values under a random assignment of observed ICs across the K runs. This ''null'' distribution is used to compute reproducibility p-values for each observed matched component from RAICAR. An objective cutoff such as pv0:05 can be used to detect ''significantly reproducible'' components. This avoids subjective user decisions such as selection of the number of clusters in ICASSO or the reproducibility cutoff in RAICAR or a cutoff on intra cluster distance in sogICA.

Results for publicly available rsfMRI data
We applied RAICAR-N to publicly available N~23 subject rsfMRI data from http://www.nitrc.org/. We analyzed the data in 2 different ways: 1. n C~4 0 ICs were extracted for each of the N~23 subjects.
The K~23 single subject ICA runs were subjected to a RAICAR-N analysis (after registration to standard space). In single subject ICA based RAICAR-N analysis (Figures 8, 9), we are able to declare 6 out of the 8 ICs reported in [8] (which used group ICA) as ''reproducible'' (4 ICs have p-values v0:05 and 2 ICs have p-values v0:06). This is consistent with the 5 reproducible RSNs reported in [32] using single subject ICA analysis.
2. L~5 subjects were randomly drawn from N~23 subjects to create one group of subjects which was subjected to a group ICA analysis in which n C~4 0 components were extracted. This process was repeated K~50 or 100 times and the resulting group ICA runs were subjected to a RAICAR-N analysis.
In group ICA based RAICAR-N analysis (Figures 10, 11, 12, and 13), we are able to declare all 8 components reported in [8] as ''reproducible'' (at pv0:05). Some of the ICs detected as ''reproducible'' in the group ICA based RAICAR-N on human rsfMRI data are not shown in [8] but do appear in the more recent paper [31]. RAICAR-N results for K~50 are almost identical to those for K~100 suggesting that K~50 runs of group ICA are sufficient for a RAICAR-N reproducibility analysis.

Single subject ICA vs Group ICA
Based on our results, it appears that single subject ICA maps are less reproducible compared to group ICA maps as illustrated in Figures 8  and 10. A single subject ICA based analysis is more resistant to subject specific artifacts. On the other hand, a group ICA based analysis makes the strong assumption that ICs are spatially identical across subjects. If this assumption is true, group ICA takes advantage of temporal concatenation to constrain the ICs spatially across subjects thereby reducing their variance. Hence, when there are no gross artifacts in individual rsfMRI data sets, group ICA is expected to be more sensitive for reproducible IC detection. As seen in Figures 9 and 11, our results agree with this proposition. All ICs declared as ''reproducible'' in the single subject based RAICAR-N analysis continue to remain ''reproducible'' in the group ICA based RAICAR-N analysis.

How should subjects be grouped for group ICA?
This raises the question of how the subjects should be grouped together for individual group ICA runs in preparation for RAICAR-N. If all N subjects are used in all group ICA analyses then there is no diversity in the individual group ICA runs. In this case, a RAICAR-N analysis will capture algorithmic variability due to non-convexity of ICA objective function but not dataset variability. Hence, our conclusions might not be generalizable to a different set of N subjects.
Another option is to randomly select L subjects out of N for each group ICA run and submit the resulting K group ICA runs to RAICAR-N. In this case, we will account for both algorithmic and data set variability via a RAICAR-N analysis. In other words, we will be able to determine those ICs that are ''reproducible'' across different sets of L subjects and across multiple ICA runs. A key question is: How should we choose L and K? In section 2.5, we proposed a simple method to determine the number of subjects L to be used in a single group ICA run out of the N subjects -the key idea is to form groups with enough ''diversity''. Multiple such group ICA runs can then be submitted to a RAICAR-N analysis for reproducibility assessment.  [8]. We are able to declare 4 ''standard'' RSNs as significantly reproducible at a p-value v0:05. There are 2 other ''standard'' RSNs that achieve a reproducibility p-value between 0:05 and 0:06 as well as 2 ''non-standard'' RSNs that achieve p-values of 0:0125 and 0:05699 respectively. We also could not find 2 of the published RSNs in [3] as reproducible in single subject ICA runs. doi:10.1371/journal.pone.0027594.g009 Clearly, the larger the value of N, the larger the value of L. Hence, increasing the number of subjects N in a study will allow us to make conclusions that are generalizable to a larger set of L subjects. Also, conclusions generalizable to L 1 subjects are expected to hold for L 2 wL 1 subjects but not vice versa.

RAICAR-N for group comparisons of reproducible ICs
In the present work, our focus was on enabling the selection of reproducible ICs for a given single group of subjects. However, RAICAR-N can be extended for between group analysis of reproducible components as well. Before we describe how to do so, it is useful to discuss other approaches for group analysis of RSNs described in Appendix S1. Suppose we have two groups of subjects A and B.
4.4.1 Discussion of single group ICA based approaches. 1. Subject specific maps corresponding to group ICA maps derived using ICA back projection or dual regression are not true ICs, i.e., they are not solutions to an ICA problem.
2. These approaches do not account for either the algorithmic or the data set variability of an ICA decomposition. The single group ICA decomposition will contain both reproducible and nonreproducible ICs, but there is no systematic way to differentiate between the two.
3. Both dual regression and ICA back projection using data derived IC templates are circular analyses. First, group ICA using all data is used to derive template IC maps or template time courses. Next least-squares based ICA back projection or dual regression using a subset of the same data is used to derive subject specific maps and time courses corresponding to each IC. Thus model 1 (group ICA) on data D is used to learn an assumption A (template IC maps or template time courses) that is then used to fit model 2 (dual regression or ICA back projection) on a subset of the same data D. This is circular analysis [33,34].
It is easy to avoid circular analysis in a dual regression approach via cross-validation. For example, one can split the groups A and B into two random parts, a ''training'' set and a ''test'' set. First, the ''training'' set can be used to derive template IC maps using group ICA. Next, the ''training'' set based template IC maps can be used as spatial regressors for dual regression on the ''test'' set. Alternatively, the template ICs for dual regression can also come from a separate ICA decomposition on a independent data set unrelated to groups A and B such as human rsfMRI data. This train/test approach cleanly avoids the circular analysis problem. It is not clear how to use cross-validation for an ICA back projection approach since template time courses cannot be assumed to remain the same across ICA decompositions.
4. Subject specific structured noise is quite variable in terms of its spatial structure. Hence, a group ICA analysis cannot easily model or account for subject specific structured noise via group level ICs. Consequently, subject specific spatial maps in ICA back projection or dual regression will have a noise component that is purely driven by the amount of structured noise in individual subjects. On the other hand, a single subject ICA based analysis can accurately model subject specific structured noise via single subject ICs.
4.4.2 Discussion of multiple ICA run approaches. 1. [35] report that using different sets of template ICs in template based methods using spatial correlation such as [36] can result in the selection of different ICs in individual ICA runs. This is not surprising since IC correspondence derived from template based methods does depend on the particular template used. This is similar to a seed based correlation analysis being dependent on the particular seed ROI used. It is worth noting that template free approaches such as sogICA and RAICAR do not rely on any template.
2. [22] state that individual runs across subjects (or groups of subjects) can be quite variable in terms of the spatial structure of the estimated ICs. For example [22], point out that an IC might be apparently split into two sub-components in some subjects but not others. The real problem is that the same model order could lead to over-fitting in some subjects (or groups of subjects) but not in others. Hence, the observed differences in a group comparison might be biased by the unknown difference in the amount of overfitting across groups A and B.    Multiple group ICA runs using groups of subjects with enough ''diversity'' can be used to account for the run-to-run variability in ICA algorithms both due to the non-convex ICA objective function as well as across subjects data variability. These group ICA runs can be subjected to a RAICAR-N ''reproducibility'' analysis. RAICAR-N enables the objective detection of ''reproducible components'' in any component based analysis of fMRI data such as ICA and can also be used for a between group comparison of ''reproducible'' ICs.

Supporting Information
Appendix S1 Background on group comparison of ICA results. (TEX)