The BFNs extraction has been formulated as a source separation problem, based on the functional integration property of the brain (McKeown et al., 1998; Du and Fan, 2013; Shi et al., 2015a). This source separation problem can usually be divided into blind source separation (BSS) and semi-blind source separation (SBSS), depending on whether the prior is given or not. With the respect to ICA model, as a representative of BSS, it is assumed that the observed fMRI mixtures (denoted as X) are linearly mixed by a set of non-Gaussian sources, namely BFNs (denoted as S), which can be formulated as

The goal of ICA is to estimate an unmixing matrix W, such that the estimated sources Y computed by the following equation are good approximation of the true sources S:

To solve Equation (1), many ICA algorithms have been proposed, e.g., the commonly used Infomax (Bell and Sejnowski, 1995) and FastICA (Hyvärinen and Oja, 1997).

ICA-R model of SBSS incorporates the prior spatial reference (denoted as r), and can be modeled in a constrained ICA framework, to the following constrained optimization problem

where J(y) is the contrast function of a standard ICA algorithm, and g(y) = ε(y, r)−ξ, with ε(y, r) denotes the closeness between y (the estimated BFN) and the reference signal r, and ξ signifies a threshold parameter used to restrain the distance between y and r. To solve Equation (3), the Lagrange multiplier method can be utilized to search for the solution using Newton-like learning (Lu and Rajapakse, 2006), fixed-point learning (Lin et al., 2007), or multi-object optimization (Du and Fan, 2013).

Supposing that there are m scanned fMRI datasets, i.e., Dataset_k, 1 ≤ k ≤ m, the fMRI data of subject i in Dataset_k is denoted as Subki, where 1 ≤ i ≤ n_k, 1 ≤ k ≤ m, and n_k signifying the subject number in Dataset_k. As described in Figure 1, the FMICA model mainly consists of three levels of ICA decomposition and two re-estimation of group-specific and subject-specific feature maps using ICA-R: (1) the first (or single-subject) level ICA decomposition on Subki to obtain the feature maps, i.e., independent components ICSki, where 1 ≤ i ≤ n_k and 1 ≤ k ≤ m; (2) the second (or intragroup) level ICA decomposition on the aggregated feature maps, i.e., ICSkAgg for Dataset_k, to obtain the feature maps at intragroup level, i.e., GICS_k for Dataset_k, 1 ≤ k ≤ m; (3) the third (or intergroup) level ICA decomposition on the aggregated feature maps across the datasets (Dataset_k, 1 ≤ k ≤ m), i.e., GICS1:mAgg, to extract the intergroup feature maps across the different datasets, i.e., GICS^; (4) the ICA-R algorithm first runs on the GICS_k regrading Dataset_k (1 ≤ k ≤ m) to extract the correspondingly intragroup-specific feature maps, i.e., GICS~k, then on the ICSki regarding Subki (1 ≤ i ≤ n_k, 1 ≤ k ≤ m) to obtain the correspondingly subject-specific featured maps (denoted as ICS~ki).

Specifically, in order to make the last procedure in FMICA more apparent, the corresponding details are described in the following. On the one hand, for extracting the intragroup-specific feature maps, i.e., GICS~k, 1 ≤ k ≤ m, the ICA-R algorithm (Du and Fan, 2013) is implemented on each GICS_k, where the correspondingly intergroup feature maps GICS^ are used as the spatial references. On the other hand, in order to extract the subject-specific feature maps (ICS~ki) corresponding to the subject data Subki, 1 ≤ i ≤ n_k, and 1 ≤ k ≤ m, the similar ICA-R procedure with spatial reference is implemented on each ICSki. It is noteworthy that the spatial references may have two choices: the intergroup feature maps GICS^ or the corresponding intragroup-specifc feature maps GICS~k. Repeating the above procedure for Dataset_k (1 ≤ k ≤ m), all the corresponding intragroup-specific GICS~k and subject-specific feature maps ICS~ki (1 ≤ i ≤ n_k) are retrieved. However, when the number of the involved datasets is <2, i.e., m = 1, the third level ICA decomposition procedure is not implemented, and both the intergroup and intragroup-specific feature maps are not generated. Facing this situation, the subject-specific feature maps (ICS~1i) are also obtained by ICA-R, where the intragoup feature maps (GICS₁) determined by the second level ICA decomposition procedure are used as the required spatial references.

Further, the corresponding statistical parametrical maps (SPMs) for ICS~ki, GICS~k, GICS^, and GICS_k (1 ≤ i ≤ n_k and 1 ≤ k ≤ m) are obtained by the z-score transformation, where the BFNs are generated by threshing the corresponding SPMs with cluster size controlling.

Finally, in FMICA, it is worthy of noting that the intragroup or intergroup BFNs are identified by the second or third level ICA decomposition procedure, while the intragroup-specific or subject-specific BFNs are both identified by ICA-R procedure using the ones at the intragroup or intergroup level as references. It is expected that the BFNs at the intergroup, intragroup-specific and subject-specific levels, have some degree of similarity to each other in spatial distribution, but orderly capture the commonness across different groups, the specific activation parts of the spatial distribution within a certain group and within a single subject data. In a word, for a given intergroup BFN, the corresponding intragroup-specific or subject-specific one belongs to the same kind of BFN, but captures the group-specific or subject-specific difference in spatial distribution of BFN.

Based on the above description, FMICA is a quite generalized framework using feature maps, which is effective to capture the common BFNs (i.e., GICS^, GICS_k), subject-specific ones (i.e., ICS~ki) and intragroup-specific ones (i.e., GICS~k). This implies that FMICA can be used to not only explore the subject-specific differences within a group, but also to reveal the intragroup-specific differences across the multiple datasets.

With respect to the FMICA implementation, the number of independent components (IC) in ICA decomposition at different levels should be first addressed. For the first level ICA decomposition procedure (depicted in Figure 1), the Laplace approximation (Minka, 2000) previously used in probabilistic ICA for fMRI data analysis (Beckmann and Smith, 2004) is used to estimate the components number for each single-subject fMRI data; for the second level ICA decomposition, the mean order corresponding to all subjects within the same dataset is used, while the average order is usually used as the number of intragroup level components in TCGICA (Calhoun et al., 2001; Li et al., 2007) implemented in the GIFT software (http://mialab.mrn.org/software/gift/index.html); for the third level ICA decomposition, the average number of components in GICS_k, 1 ≤ k ≤ m is used due to the situation of the different session scans of the same subjects under the same condition (for example, Experiment 2 of Section Experimental Designs); otherwise, the stability measure retrieved by ICASSO (Himberg et al., 2004) is used to determine the optimal components number, where ICASSO runs from the minimum number (i.e., min{x|x = #(GICS_k), 1 ≤ k ≤ m}) to the maximum one (i.e., max{x|x = #(GICS_k), 1 ≤ k ≤ m}) to obtain the stability values under different order number. Specifically, #() is an operation of obtaining the components number of the intragroup feature maps GICS_k, 1 ≤ k ≤ m (for example, Experiment 3 of Section Experimental Designs). Moreover, in this research, FMICA takes advantage of the FastICA (Hyvärinen and Oja, 1997) and GIG-ICA (Du and Fan, 2013), to perform the ICA decomposition and ICA-R decomposition, respectively. Additionally, since the performance of ICA-R depends on the accuracy of the spatial references (Du and Fan, 2013; Wang et al., 2014; Shi et al., 2015a), slightly thresholded feature maps which are more similar to the real activated BFNs are used as spatial references, where the corresponding z threshold value is set to 1.0 empirically. Finally, the z-threshold value and the cluster size threshold are set to 2.0 and 10 voxels, respectively, to obtain the ultimate BFNs.

All computations of this study were performed on a personal computer with intel(R) Core(TM) i5-3210M 2.5 GHz CPU and 4 GB RAM. The operation system platform was Windows 7. All steps for preprocessing or processing were run on the Matlab platform (Matlab 2012b, Mathworks Inc., Sherborn, MA, USA).

No preprocessing step was involved for the simulation dataset; For the other real data, the widely used DPARSF (Yan and Zang, 2010) batch processing pipeline with embedding SPM8 software (http://www.fil.ion.ucl.ac.uk/spm/) was used to perform the preprocessing operations including slice-timing, motion correction, spatial normalization to the Montreal Neurological Institute (MNI) EPI template and spatial smoothing with the full width at half maximum (FWHM) equal to 6 mm. Specifically, considering magnetization equilibrium, the first ten volumes were discarded for the test-retest NYU datasets. For the task-related datasets, no first volumes were discarded.

The z threshold of the z-scored SPMs from all the real fMRI datasets was set to 2.0, and the least active cluster size was set to 10 voxels. The BFNs were displayed by the MRIcroN software (https://www.nitrc.org/projects/mricron), and their locations were assessed by the PickAtlas toolbox (Maldjian et al., 2003, 2004).

Three kinds of experiments were designed to validate the effectiveness of FMICA in this study.

Experiment 1: the simulation dataset with only one session was used to validate effectiveness of FMICA in two aspects, i.e., the BFN identification ability at the subject-specific and intragroup levels. Specifically, the third level ICA step was not involved in this experiment. The corresponding pipeline consisted of three procedures: the first level ICA on simulation dataset for obtaining the initial FMs (ICs), the second level ICA for extracting the second level (or intragroup) FMs and the ICA-R procedure using the intragroup FMs as the references for obtaining the subject-specific BFNs, respectively.

Experiment 2: the NYU resting-state datasets with three sessions were used to validate effectiveness of FMICA in three aspects, i.e., the BFN identification ability at the subject-specific, intragroup-specific and intergroup levels. The corresponding pipeline consisted of the following steps: the first level ICA on fMRI datasets of each rest session for obtaining the initial FMs, the second level ICA for extracting the second level FMs, the third level ICA for obtaining the intergroup FMs and the ICA-R using the intergroup FMs as the references for obtaining the intragroup-specific and subject-specific FMs, respectively.

Experiment 3: the ability of capturing the group difference of intrinsic BFNs across the multiple kinds of datasets using FMICA was validated using a combination of the test-retest NYU resting-state datasets, the test-retest task-related datasets for motor, language, and spatial attention and the visual task dataset. The corresponding pipeline consisted of the following steps: the first level ICA on each aforementioned dataset for obtaining the initial FMs, the second level ICA for extracting the intragroup FMs, the third level ICA for retrieving the intergroup FMs and the ICA-R procedure using the intergroup FMs as the references for obtaining the intragroup-specific FMs.

According to Experiment 1, the BFN detection ability of FMICA at the subject-specific and intragroup levels was designed to be validated on the simulation dataset. The order in both first (individual) and second (intragroup) level ICA decomposition was set to 13 (twelve designed sources and one background source). Firstly, the 12 sources determined by FMICA at intragroup level were displayed in Figure 2, which were highly approximate to the simulated ground truth sources. The Pearson correlation coefficients between the 12 estimated intragroup sources and the corresponding 12 ground truth sources were 0.9783, 0.9672, 0.989, 0.9858, 0.9634, 0.9687, 0.9687, 0.9877, 0.9717, 0.9670, 0.9868, and 0.9703, respectively, quantitatively implying the effectiveness of FMICA in the intragroup BFNs identification. Moreover, in order to investigate the intragroup BFNs estimation from the different ratios of the retained ICA components of the intragroup level to that of the individual level, FMICA was performed with a variety of such intragroup-to-individual ratios on the simulation dataset to obtain the intragroup-level BFNs, and then the mean and standard deviation (std) values of Pearson correlation coefficients between the estimated intragroup sources and the corresponding ground truth sources at each run were calculated. As shown in Table S1, from which, one could draw a conclusion that increasing the number of retained components at the individual level had no effect on performance, while greatly increasing the number of retained components at intragroup level had a certain degree of negative impact on the estimated intragroup sources, possibly due to the over-spilt effects in ICA decomposition in the simulation dataset. Table S1 also demonstrated that using the 13 components in both first-level and second-level ICA exhibited good BFNs identification performance in this simulation dataset.

The subject-specific BFNs for each simulated subject were also estimated by FMICA, and then the correlation coefficients between these estimated subject-specific BFNs and the corresponding ground truth sources for each subject were also calculated, with the mean correlation coefficient across the 12 sources for each subject denoted as its subject-specific BFNs identification accuracy, which was compared to those of the traditional ICA model, demonstrating superior subject-specific BFNs identification ability as shown in Figure 3.

In this experiment, the test-retest NYU resting-state datasets with three sessions were used to validate effectiveness of the proposed FMICA model on identifying the BFNs at the intergroup, intragroup-specific, and subject-specific levels. At the subject-specific level, compared to the intrinsic BFNs from different subjects, the intrinsic BFNs originated from the different sessions of the same subject should be more correlated. Moreover, the intragroup-specific intrinsic BFNs from different sessions should be also highly correlated with each other due to the high reproducibility of the intrinsic BFNs (Shehzad et al., 2009; Zuo et al., 2010; Wang et al., 2016b).

Eighteen intrinsic BFNs at the intergroup and intragroup-specific levels for the resting session 1 (S1), session 2 (S2), and session 3 (S3) were selected visually by the experts from the estimated components by FMICA, as displayed in Figure 4, respectively, with the involved Talairach Daemon (TD) lobes, Brodmann areas, Automated Anatomical Labeling (AAL) atlas regions, and the representative MNI coordinates, presented in Table 1. The components IC1-IC4 were referred to the well-known default mode network (DMN; Raichle et al., 2001; Damoiseaux et al., 2006; De Luca et al., 2006), which were divided into four sub-networks (Zuo et al., 2010); IC5, IC6, IC7, and IC10 were called the auditory network, predominant visual network, lateral visual network, and sensorimotor network, respectively (Beckmann et al., 2005; Damoiseaux et al., 2006; De Luca et al., 2006; Schöpf et al., 2010; Wang et al., 2012); IC8 and IC9 were involved with the working memory function related brain regions (Wang et al., 2012, 2013; Iraji et al., 2016); IC11 and IC12 involved dorsal parietal and lateral prefrontal cortex, which were two split separate components of a dorsal pathway network (Damoiseaux et al., 2008; Schöpf et al., 2010; Wang et al., 2012, 2013); IC13 was the salience network as reported by Menon and Uddin (2010) and Uddin (2015); IC14 was the basal ganglia network, involving mainly caudate nucleus and putamen (Iraji et al., 2016); IC15 involved the cerebellum posterior lobe and a portion of calcarine area in the occipital lobe; IC16 was located at the brainstem and cerebellum, e.g., cerebellar vermis; IC17 involved mainly brodmann areas 47 and 34, e.g., the superior temporal pole; IC18 was located at a portion of the limic and frontal cortex, e.g., hippocampus, some areas of superior frontal gyrus, etc. The successful identification of the aforementioned well-known intrisic BFNs at the intergroup and intragroup-specific levels demonstrated the effectiveness of the proposed FMICA model.

From Figure 4, it could be observed that the intrinsic BFNs at the intragroup-specific level from each session were quite approximate to the corresponding ones at the intergroup level. Meanwhile, three pairs of mutual correlations among the intragroup-specific BFNs estimated from the three sessions, were calculated, as shown in Figure 5, from which high correlations could be clearly observed, demonstrating great reproducibility of the intrinsic BFNs across the sessions.

At the subject-specific level, the BFNs identified by FMICA were compared among sessions of the same subject and among different subjects, respectively, and mean correlation coefficients between pairs of BFNs under comparison for each subject were shown in Figure 6, where Figures 6A,B were for across-sessions and across-subjects comparison, respectively. Similar to FMICA, the BFNs identified by the first level ICA were also compared among sessions of the same subject and among different subjects, respectively, and mean correlation coefficients between pairs of BFNs under comparison for each subject were also shown in Figures 6C,D for across-sessions and across-subjects comparison, respectively. It was worth noting that the 18 intrinsic BFNs at the intergroup level were used as the templates to match the best ones from the individual FMs generated by the first level ICA decomposition on each session data of each subject, aiming at overcoming the random order of the components. Moreover, based on the contrast values presented in Figures 6A,C, and the contrast ones in Figures 6B,D, two sample T-tests with significance level equal to 0.05 were implemented, respectively, where the mean value (0.5644, marked in Figure 6A) of all points in Figure 6A was significantly larger than the one (0.3168, marked in Figure 6C) of all points in Figure 6C with p value equal to 6.1630 × e⁻⁸⁰, and the mean value (0.3591, in Figure 6B) of all points in Figure 6B was also significantly larger than the one (0.1853, marked in Figure 6D) of all points in Figure 6D with p value equal to 2.1511 × e⁻¹⁷². The mean values of the points in Figures 6B,D corresponding to the first level ICA were relatively small, and this was possibly due to that some BFNs could be identified at the intergroup level, but not separated at the single subject level by the traditional ICA. However, the ICA-R re-estimation procedure in FMICA could identify most of BFNs at the single subject level, demonstrating that the proposed FMICA could identify the subject-specific BFNs more effectively than the traditional ICA did.

To intuitively compare the performance between the proposed FMICA and the traditional ICA, the spatial maps of IC1 (DMN) identified by FMICA and FastICA for all three sessions of the first four subjects (for space limitation) as shown in Figures 7A,B, respectively, were taken as an example, from which it was obvious that the DMNs identified by FMICA had higher across-sessions than across-subjects consistency and much higher both across-sessions and across-subjects consistency compared to the results of FastICA, implying that FMICA had higher subject-specific BFNs identification capability than traditional ICA.

In summary, results from the test-retest resting-state datasets demonstrated that the proposed FMICA model had high BFN identification capability at intergroup, intragroup-specific and subject-specific levels.

In this experiment, the intergroup and intragroup-specific analysis ability of FMICA was further validated by combining the test-retest resting-state datasets, the test-retest task-related datasets for motor, language and spatial attention (i.e., Task1, Task2, and Task3) and the visual task dataset. Namely, there were ten datasets as the input of FMICA, i.e., the resting-state datasets with three sessions, three test-retest task-related datasets (i.e., Task1, Task2, and Task3) with two sessions and the visual task dataset with one session. In this experiment, as described in Section Some Key Points in FMICA Implementation, the ICASSO method was used to determine the optimal order for the intergroup-level analysis based on the stability measure of the estimated components for each order, as shown in Figure S2, demonstrating that the estimated components had the highest mean/median stability and relatively small values of standard deviation (STD) and inter-quartile range (IQR), when the order was equal to 57. Therefore, the order was set to 57 in Experiment 3.

Intragroup-specific BFNs for each of the 10 datasets and the corresponding intergroup BFNs, selected visually by the experts from the estimated components by FMICA, were showed in Figure 8 for the first five BFNs due to space limitation, and results for the remaining 25 BFNs were shown in Figure S3. The TD lobes, Brodmann areas, AAL regions and the representative MNI coordinates involved in these BFNs, were recorded in Table S2. It could be observed that most of the intrinsic BFNs extracted in Experiment 2 had high reproducibility in Experiment 3, and better across-sessions than across-datasets consistency of the BFNs was also observed. Correlation analysis was performed to quantify such consistency in various cases. Firstly, mean and std values of correlation coefficients among the intragroup-specific BFNs from different sessions of the resting-state datasets under the same condition (e.g., session 1 and session 2 of resting-state datasets), were calculated and presented in Figure 9A, demonstrating a high mean correlation of 0.8464 and thus implying high across-sessions reproducibility of the BFNs in resting-state datasets (Wang et al., 2016b). Then, the same correlation analysis were performed across the test-retest task datasets of Tasks 1, 2, and 3 from the same subjects, with results presented in Figure 9B, showing also a high mean correlation of 0.8314 and thus implying across-tasks similarity of the intrinsic functional connectivity architecture (Finn et al., 2015). Finally, correlation analysis were performed on completely different kinds of datasets sharing neither sessions nor tasks, with results shown in Figure 9C, indicating a low correlation of 0.5814 inferior to that in Figures 9A,B and thus demonstrating that the proposed FMICA was able to effectively capture the differences of intragroup BFNs across the different kinds of datasets.

To summarize, it could be stated that the proposed FMICA was effective for the intergroup and intragroup-specific analysis, and could characterize the group-specific difference.

It was proposed to perform ICA analysis by summarizing fMRI data of each subject as a feature map and applying subsequently traditional ICA algorithms on these feature maps, where features could be the amplitude of low frequency fluctuations (ALFF) maps for resting-state data or T-statistic maps for task-related data, yielding BFNs strikingly similar to but slightly noisier than the results of spatiotemporal group ICA analysis (i.e., TCGICA; Calhoun and Allen, 2013). Very recently, another novel feature-based ICA model using seed-based functional connectivity as summarizing features was proposed (Iraji et al., 2016), with performance highly depending on the choice of seeds. The proposed FMICA in this paper took spatial maps of the independent components of the subjects as feature maps for group analysis. FMICA could produce the comparable intragroup BFNs to the spatiotemporal-domain based group ICA, as shown in Figure S4, implying that it was more effective than the first feature-based ICA model. Meanwhile, FMICA without the seed-based functional connectivity identification procedure was more flexible than the second feature-based ICA model, and it had extra unique advantages of identifying subject-specific, intragroup-specific and intergroup BFNs due to hierarchical processing incorporated in the model.

Single-subject independent components were used as the input feature maps in the proposed FMICA model. However, edges and shapes of the feature-maps could be susceptible to the preprocessing steps in fMRI data analysis, such as the spatial smoothing with FWHM kernels of different size. Therefore, one future research topic on FMICA might be to develop a more robust model to deal with the effects of the preprocessing steps on the feature maps.

Since brain activity at either resting or task states is non-stationary, and it is very important to characterize the dynamics of brain networks (Calhoun et al., 2014). Although static brain functional activity is considered in the current study, the FMICA has also the potential to provide new options to the investigation of the dynamic characteristics of brain networks.

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