Group-ICA Model Order Highlights Patterns of Functional Brain Connectivity.

Resting-state networks (RSNs) can be reliably and reproducibly detected using independent component analysis (ICA) at both individual subject and group levels. Altering ICA dimensionality (model order) estimation can have a significant impact on the spatial characteristics of the RSNs as well as their parcellation into sub-networks. Recent evidence from several neuroimaging studies suggests that the human brain has a modular hierarchical organization which resembles the hierarchy depicted by different ICA model orders. We hypothesized that functional connectivity between-group differences measured with ICA might be affected by model order selection. We investigated differences in functional connectivity using so-called dual regression as a function of ICA model order in a group of unmedicated seasonal affective disorder (SAD) patients compared to normal healthy controls. The results showed that the detected disease-related differences in functional connectivity alter as a function of ICA model order. The volume of between-group differences altered significantly as a function of ICA model order reaching maximum at model order 70 (which seems to be an optimal point that conveys the largest between-group difference) then stabilized afterwards. Our results show that fine-grained RSNs enable better detection of detailed disease-related functional connectivity changes. However, high model orders show an increased risk of false positives that needs to be overcome. Our findings suggest that multilevel ICA exploration of functional connectivity enables optimization of sensitivity to brain disorders.

Introduction

Magnetic resonance imaging (MRI) has developed rapidly during recent years, enabling very accurate structural and functional inferences of neurological and psychiatric diseases. Functional MRI (fMRI) enables the detection of task responses as well as spontaneous interregional connectivity assessment of the living human brain without invasive or radioactive methodology. Resting-state functional connectivity analyses study similarities in the temporal behavior of blood oxygen level dependent (BOLD) signal fluctuations in different brain regions (Biswal et al., 1995; Cordes et al., 2000; Lowe et al., 2000; Greicius et al., 2003; Beckmann et al., 2005; Fox et al., 2005). Coherent spatial patterns of low-frequency (<0.1 Hz) fluctuations in the resting-state BOLD signal are referred to as a functional network. Alterations in functional connectivity of such networks are suggested to precede both structural changes and clinical symptoms (Greicius et al., 2004; Filippini et al., 2009). A number of studies have used independent component analysis (ICA) approaches to measure functional connectivity in clinical populations such as Alzheimer's disease or dementia (Greicius et al., 2004; Rombouts et al., 2009; Seeley et al., 2009), schizophrenia (Jafri et al., 2008; Calhoun et al., 2009), depression (Anand et al., 2005; Greicius et al., 2007; Chen et al., 2008; Zhou et al., 2009; Sheline et al., 2010), epilepsy (Zhang et al., 2009), Huntington's disease (Wolf et al., 2008), and amyotrophic lateral sclerosis (Mohammadi et al., 2009).

Independent component analysis as a blind source separation technique has become a major data-driven analysis tool for fMRI studies (McKeown et al., 1998; Biswal and Ulmer, 1999; Calhoun et al., 2001; Kiviniemi et al., 2003). It is an explorative data analysis method that produces a number of spatial maps (spatial components) and corresponding time courses (Calhoun et al., 2001). Generally, spatial ICA is a more appropriate method than time-domain ICA for analyzing resting-state fMRI data, given the small number of time points and large number of voxels (spatial samples) included in most fMRI datasets. Resting-state networks (RSNs) can be reliably and reproducibly detected using ICA at individual subject and group levels (Greicius et al., 2004; Damoiseaux et al., 2006; Shehzad et al., 2009; Zuo et al., 2010). Generally, the identification of meaningful neurophysiological spatial components is usually performed either by spatial correlation with a predefined template (Greicius et al., 2003; Van de Ven et al., 2004; Calhoun et al., 2008) or by visual inspection (Damoiseaux et al., 2006; DeLuca et al., 2006; Harrison et al., 2008).

Importantly, altering the dimensionality (model order) estimation in ICA can have a significant impact on the spatial characteristics of the RSNs identified (Abou Elseoud et al., 2010). Accordingly, ICA results may be “split” into a number of sub-networks, depending on the parameters of the analysis (e.g., model order selection). Notably, the process of ICA model order selection is somewhat arbitrary (i.e., one has to tell ICA how many components to estimate), depending on a number of factors, e.g., data quality (Strother et al., 2002, 2010), time points i.e., if one does not perform data reduction, the interpretation of the component maps can be a very time consuming task and a large majority of these maps are not useful (Calhoun et al., 2004), field strength (as the signal-to-noise ratio of BOLD signal increases with field strength), number of subjects (Yourganov et al., 2010), and regions or functions of interest (Kiviniemi et al., 2009; Abou Elseoud et al., 2010). In addition, automatic model order estimation may not be reliable enough to be implemented as a standard methodology (Yourganov et al., 2010), especially when comparing different studies of a given disease.

The decomposability of a network can be measured using modularity (Guimerà et al., 2004; Newman and Girvan, 2004), which can be used as a merit function to find the optimal partition of a network. There is strong evidence of brain modularity (Bullmore and Sporns 2009), arising from recent human neuroimaging studies showing anatomical (Chen et al., 2008) and functional (Ferrarini et al., 2009; Meunier et al., 2009a, 2009b) evidences for modularity of brain networks. Animal studies have supported such hierarchical organization (Hilgetag et al., 2000; Schwarz et al., 2008). The underlying functional mechanisms of the modularity of brain networks could be explained by the free-energy principle (Friston, 2009). Brain modularity, as shown by the free-energy principle, is essential for transmitting prediction errors to higher cortical areas, which use these errors to update an internal model that generates top–down predictions of sensory inputs (Friston, 2009, 2010). Meunier et al. (2009b) applied a computational algorithm to derive a hierarchical modular decomposition of human brain networks using fMRI. Eight large modules were depicted at the highest level of the hierarchy, each compromising more than 10 nodes. While at the lowest level, there were 57 sub-modules. Notably, these results stand very similar to ICA decompositions obtained at low model orders (where large-scale networks represent the large modules) as well as at high model orders (where fine-grained sub-networks represent the sub-modules).

Independent component analysis studies in clinical populations have reported functional connectivity differences, but the use of different ICA model orders makes the comparison of these results difficult. Previously, we have shown that how RSNs’ characteristics change as a function of ICA model order concerning a population of healthy subjects. There are strongly independent components (e.g., secondary sensory motor and basal ganglia) that cannot be depicted at low model orders, and, leaving them out from the analysis might lead to false negative results. Moreover, at higher model orders IC sources are finely clustered and therefore might be more sensitive to subtle connectivity alterations.

Based on our previous findings, we hypothesize that the detected disease-related differences in functional connectivity alter as a function of ICA model order. In order to investigate this hypothesis, we investigated between-group differences in functional connectivity with dual regression technique (Beckmann et al., 2009; Filippini et al., 2009). The null hypothesis was that the total volume of the detected between-group differences would not be affected by ICA model order and presents a straight linear relationship. We utilized seasonal affective disorder (SAD) as an example of a neuropsychiatric disease to compare with normal healthy controls (HCs). Our findings in SAD show only increased connectivity in SAD, which facilitates inferences on the effects of model order selection. Finally, functional connectivity changes in SAD involve RSNs that could be easily identified and followed up throughout estimated model orders.

Materials and Methods

Participants

This research is part of the SAD and light therapy project that has started in 2009 and still going on at Oulu University, Finland. In our study, SAD patients were recruited through advertisements in two waves during January–February 2009 (first wave) and November 2009–January 2010 (second wave) in the city of Oulu, Finland (latitude 65°01′N). The first wave represents the pilot study of the SAD and light therapy project, while the second wave represents the continuation of the same project. All SAD patients (39.78 ± 10.64 years, 30 ♀, 15 ♂) were interviewed by an experienced psychiatrist. Diagnostic and Statistical Manual of Mental Disorders (American Psychiatric Association, 1994) diagnoses for recurrent major depression (moderate or severe) were obtained using the Mini International Neuropsychiatric Interview (MINI; Sheehan et al., 1998). In addition, patients had to fulfill the diagnostic criteria for “seasonal pattern” according to DSM-IV-TR [American Psychiatric Association, 2000; although the diagnostic criteria for “seasonal pattern” can be applied to a diagnosis of major depressive episodes in both bipolar (I and II) disorder and recurrent MDD (American Psychiatric Association, 2000), only patients with recurrent unipolar depression were included in the present study to increase the homogeneity]. The ethical committee of Oulu University Hospital has approved the study for which the subjects have been recruited, and informed consent has been obtained from each subject individually according to the Helsinki declaration.

The exclusion criteria were as follows: lifetime psychotic disorder, other concurrent DSM-IV axis I except anxiety disorder, clinically significant DSM-IV axis II disorder, substance abuse or dependence, tobacco smoking, lifetime suicide attempt or suicide ideations during the past month, unstable physical disorder, psychotropic medications or corresponding herbal preparations, bright light therapy for the current episode, ocular-disorders except myopia/hyperopia. Furthermore, pregnant candidates were excluded. In addition, normal exclusion criteria for MRI-scanning were used. All SAD patients except four had no comorbid physical disorders (in the second wave: the first patient with arterial hypertension controlled by an angiotensin II receptor antagonist: the second patient suffers from arterial hypertension and hypercholesterolemia controlled by angiotensin II receptor antagonist and statins, respectively: the third patient is diagnosed with androgenic alopecia which is controlled by finasteride, and the fourth patient suffers from menopausal syndrome, using estradiol).

Imaging methods

Altogether, 45 anti-depressant-free SAD patients and 45 age-, gender- and ethnicity-matched HCs (no concomitant medications) from the general population were imaged using the same protocol during the same winter-period. Resting-state BOLD data were collected on a GE Signa 1.5 Tesla whole body system with an eight channel receive coil, using an EPI GRE sequence (TR 1800 ms, TE 40 ms, 280 time points, 28 oblique axial slices, slice thickness 4 mm, inter-slice space 0.4, covering the whole brain, FOV 25.6 cm × 25.6 cm, with 64 × 64 matrix, parallel imaging factor 2, and a flip angle of 90°). T1-weighted scans were imaged using 3D FSPGR BRAVO sequence (TR 12.1 ms, TE 5.2 ms, slice thickness 1.0 mm, FOV 24.0 cm, matrix 256 × 256, and flip angle 20°, and NEX 1) in order to obtain anatomical images for co-registration of the fMRI data to standard space coordinates. The subjects were instructed to simply lay still inside the scanner with their eyes closed, think of nothing particular and not to fall asleep. Motion was minimized using soft pads fitted over the ears and hearing was protected.

Data pre-processing

Head motion in the fMRI data was corrected using multi-resolution rigid body co-registration of volumes, as implemented in FSL 3.3 MCFLIRT software (Jenkinson et al., 2002). The default settings used were: middle volume as reference, a three-stage search (8 mm rough + 4 mm, initialized with 8 mm results + 4 mm fine grain, initialized with the previous 4 mm step results) with final tri-linear interpolation of voxel values, and normalized spatial correlation as the optimization cost function. Brain extraction was carried out for motion corrected BOLD volumes with optimization of the deforming smooth surface model, as implemented in FSL 3.3 BET software (Smith 2002) using threshold parameters f = 0.5 and g = 0; and for 3D FSPGR volumes, using parameters f = 0.25 and g = 0. This procedure was verified with visual inspection of the extraction result. In some cases when the eye/tonsil tissue was not removed appropriately, these tissues were extracted manually. The resulting image was used as a mask for a secondary brain extraction. After successful brain extraction the BOLD volumes were spatially smoothed with Gaussian kernel (7.5 mm FWHM) and voxel time series were detrended using a Gaussian linear high-pass filter with a 100 s cutoff. The FSL 4.1.4 fslmaths tool was used for these steps. Multi-resolution affine co-registration as implemented in the FSL 4.1.4 FLIRT software (Jenkinson et al., 2002) was used to co-register mean non-smoothed fMRI volumes to 3D FSPGR volumes of corresponding subjects, and 3D FSPGR volumes to the Montreal Neurological Institute (MNI) standard structural space template (MNI152_T1_2mm_brain template included in FSL). Tri-linear interpolation was used, a correlation ratio was used as the optimization cost function, and regarding the rotation parameters a search was done in the full [−π π] range. The resulting transformations and the tri-linear interpolation were used to spatially standardize smoothed and filtered BOLD volumes to the MNI standard space. However, for computational reasons pertaining to later analysis steps, 4 mm resolution was retained after spatial normalization.

ICA analysis

We have used spatial ICA in this paper and for simplicity, in the remainder of this paper, we refer to spatial ICA as ICA. ICA analysis was carried out using FSL 4.1.4 MELODIC software implementing probabilistic independent component analysis (PICA; Beckmann and Smith, 2004). Multisession temporal concatenation tool in MELODIC (implementing FastICA algorithm) was used to perform PICA related pre-processing and data conditioning in group analysis setting. ICA using 20 (low model order), 40, 60, 70, 80, 100, 120, and 150 (high model orders) independent component maps (IC maps) was applied to detect RSNs as described earlier (Kiviniemi et al., 2009; Abou Elseoud et al., 2010). The IC maps were thresholded using an alternative hypothesis test based on fitting a Gaussian/gamma mixture model to the distribution of voxel intensities within spatial maps (Beckmann et al., 2005) and controlling the local false-discovery rate at p < 0.5. Repeatability measures, e.g., ICASSO (Himberg et al., 2004), were not used in this study as our recent results (Remes et al., 2010) suggest that they provide only little improvement to IC estimates using FastICA. The expectation for the lack of differences is even more pronounced in a group-ICA setting due to increased SNR from pooling individual fMRI datasets into a single analysis.

The between-subject analysis of the resting data was carried out using a regression technique (dual regression) that allows for voxel-wise comparisons of resting-state fMRI (Beckmann et al., 2009; Filippini et al., 2009; Littow et al., 2010; Veer et al., 2010). Dual regression approach identifies subject-specific temporal dynamics and associated spatial maps within each subject's fMRI data set. This involves (A) multiple linear regression of the z-score group-PICA maps against the preprocessed individual 4D resampled data sets yielding a subject-specific variance normalized (des norm = 1) time course for each component separately, and (B) multiple linear regression of these time courses was carried out against the preprocessed individual data sets in order to obtain subject-specific spatial maps.

Statistical difference was assessed non-parametrically using permutation testing implemented in FSL's Randomize tool, Version 2.1, incorporating also threshold-free cluster enhancement (TFCE; Smith and Nichols, 2009). This involved deriving null distributions of t-values for the contrasts reflecting the between-group effects by performing 500 random permutations of group labels and testing the difference between groups for each permutation (Nichols and Holmes, 2002). For each RSN, the resulting statistical maps were thresholded at p < 0.05 (TFCE corrected for family wise errors). The resulting between-group difference maps were resampled into 2 mm.

Initially at low model order (most often used in the literature), RSNs were identified as anatomically and functionally classical RSNs upon visual inspection, ICA prominent low-frequency power on Fast Fourier Transformation (FFT) spectra, and slow fluctuation in time courses. RSNs were assigned to commonly described RSNs as previously reported (e.g., DeLuca et al. 2006; Kiviniemi et al., 2009; Smith et al., 2009; Abou Elseoud et al., 2010). At high model orders, RSNs were identified via spatial correlation coefficients (fslcc tool in FSL) using low model order RSNs as templates, and then verified by visual inspection across estimated model orders. The Juelich histological atlas (Eickhoff et al., 2007) and the Harvard–Oxford cortical and subcortical atlases (Harvard Center for Morphometric Analysis) provided with the FSL software were used to identify the anatomical characteristics of both RSNs and between-group differences.

For demonstration of the impact of ICA model order on between-group differences in functional connectivity, we have selected the motor and the visual large-scale networks. These RSNs were selected as they (a) showed most distinct disease-related changes and (b) could be further divided into several sub-networks at high model orders. The FSL fslstats tool was used to calculate the volume of non-zero voxels in each RSN and each significant between-group difference. These values were then divided by the number of voxels of the anatomical templates included in FSL in order to provide anatomical localization of the detected significant between-group differences as well as the RSNs. Graphs showing the total volume (number of non-zero voxels) of both between-group differences and RSNs as a function of ICA model order were made using Origin software (OriginPro 8 SR0, V8.0725). These total volumes were calculated by summing up all non-zero brain voxels using the FSL 4.1.4 fslmaths tool. The summed maps were selected from group-PICA maps (p < 0.5 threshold) and from between-group difference maps (TFCE corrected p < 0.05). RSN maps in Figure 2 are conservatively thresholded (z-score > 5) in order to show differences in functional segmentation across model orders. These maps show the core and the number of detected RSNs at each model order.

Figure 2

Functional segmentation of resting-state networks (RSNs) at different functional hierarchical levels superimposed on an MNI template. RSNs are thresholded at z-score > 5. Model order 20 yielded 11 large-scale networks (top). 47 and 70 fine-clustered RSNs obtained from model order 70 (middle) and 100 (bottom), respectively. The same color templates (Fslview color templates) were used to mark the fine-clustered and large-scale RSNs. Numbers at the bottom of the images refer to MNI coordinates (xyz).

Results

The results show significant increased functional connectivity in SAD compared to HCs at all estimated model orders. A straight linear fit without any change as a function of model order (that corresponds to our null hypothesis) does not fit the data (R2 = 0.009, p = 0.8) and it has to be discarded. Contrary to the null hypothesis, the total volume of between-group differences significantly increased (R2 = 0.6, p = 0.0006) according to the third polynomial fit, peaking at model order 70 and then decreasing gradually as a function of ICA model order (Figure 1B). On the other hand, the total RSNs volume gradually increased, and then showed relative stability at higher model orders (Figure 1A). Notably, the total volume of between-group differences showed a marked elevation between model orders 60 and 70.

Figure 1

(A) RSNs total volume significantly increases as function of ICA model order reaching model order 60. Notably, higher model orders 60–150 showed relatively stable total volume. In gray, the total number of brain voxels used in the ICA analysis. (B) Total volume of between-group differences (see Figures 3, 4, green color) shows a non-linear increase in volume as a function of model order with volume maxima at model order of 70.

Figure 3

SAD significant increased functional connectivity (green) in the motor cortex at different ICA model orders (20, 40, 60, 70, 80, 100, 120, and 150) shown superimposed on MNI template. One RSN (red–yellow) involving both the motor cortices and the auditory cortex as well at the low model order of 20. At higher model orders this large-scale RSN splits into smaller fine-grained sub-networks. Notably, some of these networks still show significant increased connectivity, while others do not. Numbers at the bottom of the images refer to MNI coordinates (xyz). The left hemisphere corresponds to the right side and t-score threshold is shown at the right.

Figure 4

SAD significant increased functional connectivity (green) in the visual cortex at different ICA model orders (20, 40, 60, 70, 80, 100, 120, and 150) shown superimposed on MNI template. While at low model order of 20 the visual cortex was segmented into one RSN (red–yellow) which also showed a significant increased connectivity, higher model orders provided a more fine-grained segmentation. A prominent increase in the number of sub-networks (red–yellow) with significant increased connectivity is evident at higher model orders. Numbers at the bottom of the images refer to MNI coordinates (xyz). The left hemisphere corresponds to the right side and t-score threshold is shown at the right.

Moreover, at model orders 70 and 100, the total volume of between-group connectivity differences in proportion to the total RSNs volume (fraction of significant voxels) is nearly the same ≈ 0.075 (see Figure S1 in Supplementary Material). The fraction of significant voxels is highly correspondent to the volume of between-group differences due to narrow differences between total RSN volumes across model orders (Figure 1A). In this paper, for simplicity we categorized the estimated ICA model orders into two groups: large-scale and fine-grained levels.

Large-scale RSNS (high hierarchy)

Large-scale RSNs were localized using a low ICA model order of 20. At this high hierarchical level, the entire brain was segmented into 11 large-scale RSNs (Figure 2). SAD patients showed significant increases in functional connectivity in four out of the 11 identified RSNs involving the precuneus, the visual cortex, the motor and the somatosensory cortices as well as the bilateral caudate and thalamus nuclei (Table 1). Most of the large-scale RSNs appear to involve more than one functional node, e.g., the motor RSN involves both the entire motor, premotor and somatosensory cortices, in addition to the auditory cortex.

Table 1

Large-scale resting-state networks (RSNs) which involved significant increased connectivity at model order 20.

RSN	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
1	Superior parietal lobule	56.0	7378	50	30	62	47.34	37.37	62.81
	Lateral occipital cortex	44.7	12130	50	28	62	45.72	27.98	56.26
	Precuneous cortex	44.4	8038	50	28	62	44.67	32.65	60.08
	Supramarginal gyrus	33.5	3581	64	42	56	58.69	44.18	57.44
	Lateral occipital cortex	29.7	5142	24	24	46	47.99	27.82	38.19
	Cuneal cortex	27.1	2139	40	27	62	43.42	26.65	54.26
	Inferior temporal gyrus	26.2	1994	72	30	34	45.3	33.21	30.73
	Middle temporal gyrus	25.8	2507	72	30	34	46.48	32.01	34.4
	Supramarginal gyrus	25.2	3744	64	42	56	51.08	41.38	58.56
	Postcentral gyrus	24.2	5736	50	34	65	48.2	41.88	62.36
	Angular gyrus	23.4	3201	38	33	64	45.85	37	57.86
	Cingulate gyrus, posterior division	18.1	1848	46	38	62	44.52	40.22	56.02
	Supracalcarine cortex	16.9	1059	39	26	60	42.34	30.57	46.54
	Temporal occipital fusiform cortex	14.9	1113	70	32	30	52.97	38.28	27.19
	Parietal operculum cortex	14.0	969	74	48	54	61.11	46.78	52.99
	Middle frontal gyrus	13.9	2531	58	66	64	47.87	66.66	61.67
	Superior frontal gyrus	13.1	2486	58	66	64	47.2	65.85	62.52
	Occipital pole	11.6	1804	52	23	60	50.37	19.78	50.54
	Intracalcarine cortex	10.1	721	36	32	44	43.25	32.47	41.99
	Precentral gyrus	9.2	2681	43	37	63	45.55	51.61	63.8
	Parahippocampal gyrus	9.1	374	60	42	26	56.24	44.14	25.99
2	Visual cortex V2 BA18	75.2	15141	36	18	48	44.35	21.92	39.58
	Visual cortex V1 BA17	73.5	12036	44	20	44	44.38	21.51	39.87
	Visual cortex V5 R	53.2	923	22	24	36	22.55	27.41	35.36
	Visual cortex V5 L	44.9	713	66	21	38	65.19	24.51	37.45
3	1ry auditory cortex TE1.0 R	90.7	1216	18	52	42	19.48	55.43	39.92
	1ry motor cortex BA4p R	84.3	4391	30	50	68	28.66	50.95	62.32
	1ry auditory cortex TE1.2 R	84.0	807	14	60	43	16.56	60.1	38.74
	1ry somatosensory cortex BA3b R	83.2	5519	30	50	68	25.66	52.02	60.58
	1ry somatosensory cortex BA3a R	81.8	2881	28	48	66	27.88	51.44	58.92
	1ry motor cortex BA4a R	78.8	5489	30	50	68	31.38	52.23	63.95
	2ry somatosensory cortex OP4 R	76.1	2401	16	58	51	15.75	59.19	44.26
	1ry motor cortex BA4a L	75.7	5791	46	58	62	58.39	51.63	63.8
	1ry somatosensory cortex BA1 R	74.5	4223	30	50	68	23.84	51.27	61.65
	1ry somatosensory cortex BA3b L	73.9	5668	64	50	66	63.84	51.23	60.48
	1ry motor cortex BA4p L	73.9	4376	64	50	66	60.88	50.29	62.07
	2ry somatosensory cortex OP4 L	73.4	2233	72	56	54	72.86	56.7	44.51
	1ry somatosensory cortex BA2 R	72.3	3558	26	48	65	24.5	48.64	61.19
	1ry somatosensory cortex BA2 L	72.1	4174	64	50	66	65.31	48.5	61.6
	1ry somatosensory cortex BA1 L	72.0	4776	64	50	66	65.57	50.17	61.89
	1ry auditory cortex TE1.0 L	71.3	1235	68	50	42	69.55	53.34	40.83
	1ry auditory cortex TE1.2 L	69.2	900	74	54	40	72.52	57.01	39.95
	1ry somatosensory cortex BA3a L	69.0	2855	64	48	65	61.36	50.47	58.65
	1ry auditory cortex TE1.1 R	68.5	795	18	52	42	22.12	52.11	40.98
	2ry somatosensory cortex OP1 R	66.1	1884	14	60	46	17.4	54.41	45.46
	2ry somatosensory cortex OP3 R	62.3	1420	22	56	52	20.62	57.79	43.87
	1ry auditory cortex TE1.1 L	54.6	785	68	50	42	66.84	50.8	41.09
	Premotor cortex BA6 R	52.8	6931	30	50	68	33.44	55.85	64.11
	2ry somatosensory cortex OP3 L	51.9	823	74	60	45	69.32	56.53	42.14
	2ry somatosensory cortex OP2 R	51.6	637	22	51	42	24.14	54.84	42.57
	2ry somatosensory cortex OP1 L	49.8	1833	70	55	54	71.53	52.92	45.17
	Premotor cortex BA6 L	48.7	6009	46	58	62	55.8	55.38	64.41
	Anterior intra-parietal sulcus hIP2 R	35.7	1004	24	48	64	23.41	46.7	60.84
	2ry somatosensory cortex OP2 L	35.4	296	68	50	43	65.22	52.24	43.22
	Anterior intra-parietal sulcus hIP2 L	26.8	777	68	52	58	68.58	46.39	59.28
	BA44 R	17.2	716	14	62	46	14.85	64.83	44.54
	BA44 L	13.9	934	75	58	50	74.3	61.16	44.29
	Anterior intra-parietal sulcus hIP1 R	12.3	475	26	44	64	25.67	42.64	61.53
4	Left Putamen	34.9	820	61	69	36	55.87	68.25	35.65
	BA45 L	29.4	1497	68	70	34	67.84	75.03	38.76
	Left Pallidum	28.3	309	52	66	38	52.96	64.46	36.49
	Left Caudate	27.9	528	52	66	38	51.67	68.12	37.91
	Right Putamen	27.8	695	28	70	37	33.92	69.42	35.2
	Premotor cortex BA6 R	27.0	3551	44	72	56	39.16	68.02	62.5
	Premotor cortex BA6 L	26.5	3274	46	70	58	51.39	67.31	63.56
	Left Accumbens	25.6	137	51	68	36	51.47	69.67	34.22
	Left Thalamus	23.0	653	50	65	38	48.54	56.38	37.91
	Right Caudate	22.9	490	38	68	38	37.83	68.46	37.36
	Right Pallidum	21.5	268	38	68	38	37.01	65.19	36.33
	BA45 R	18.8	812	18	71	32	20.44	73.4	37.19
	BA44 L	18.1	1215	62	82	47	68.24	72.01	39.13
	Right Thalamus	16.7	503	38	65	40	40.76	57.71	38.37
	BA44 R	16.7	693	22	70	34	20.82	68.95	41.94

Large-scale RSNs are demonstrated by: anatomical regions involved, number of voxels, and MNI coordinates (in mm) of maximum z-scores, and center of mass.

Fine-grained sub-networks (low hierarchies)

At higher model orders, large-scale RSNs branched into sub-networks segmenting every network into smaller fine-grained sub-networks. For instance, the large-scale single sensorimotor RSN was segmented into 5–11 sub-networks (at model order 40–150) covering the entire motor, premotor and somatosensory cortices, in addition to, a network involving the auditory cortex (Figure 3). The remaining large-scale RSNs showed similar tendency to split down into either right and left, anterior and posterior, or superior and inferior compartments. Segmentation of the brain functionality into detailed sub-networks using ICA model orders of 40, 60, 70, 80, 100, 120, and 150 yielded 27, 40, 47, 54, 70, 81, and 95 RSNs, respectively (Table 2). At each model order, almost every network is decomposed of a number of smaller units (see Tables 1–8 in Supplementary Material). Therefore, each of these model orders represents a unique functional hierarchical level.

Table 2

The table shows the total number of resting-state networks (RSNs) and RSNs with a significant abnormal connectivity.

No. of model order	Total no. of RSNs	No. of SAD increased connectivity RSNs
20	11	4
40	27	16
60	46	22
70	47	25
80	55	36
100	70	40
120	81	44
150	95	46

Both numbers show a linearly increase as a function of ICA model order.

Figure 2 shows functional network segmentations at both large-scale (20 ICs) and fine-grained (70 and 100 ICs) levels. RSNs are thresholded at z-score > 5 for demonstration purposes. Notably, comparing these hierarchical levels, fine-grained level RSNs are segmented into sub-network clusters, i.e., right and left, anterior and posterior, or superior and inferior, etc. (Figure 2, the motor cortex in light blue color). Indeed functional segmentation at even higher model orders, e.g., 150, is possible and will yield focused detailed clusters, but at the same time is not feasible for demonstration purposes because such large number of clusters overlap widely (therefore, model orders 20, 70, and 100 were chosen for demonstration). Importantly, complex subcortical structures, e.g., basal ganglia, which are parts of large-scale RSNs at low model orders (not shown in Figure 2, z-score > 5) are clearly depicted as separate networks at high model orders (see Figure 2, e.g., thalamus and caudate in hot color).

The effect of functional hierarchy on between-group significant differences

The results showed a significant linear increase in the total number of RSNs as a function of ICA model order. Interestingly, at low model order of 20, only one motor RSN covers the motor brain areas and also one visual RSN covers the visual areas (see Figures 3, 4). At model order 70, the total number of the motor RSNs was eight, while the visual RSNs consisted of 10 networks. At model order 150, the motor and visual RSNs consisted of 11 and 16 networks, respectively (see Figures 3, 4). Notably, some fine-grained RSNs do not involve between-group differences, although their lower model order RSN precursors show significant between-group differences, i.e., center motor and secondary somatosensory network (see Figures 3, 4 and Tables 3–7). It is obvious that the detected between-group differences are distributed differently at each model order, particularly at the highest estimated model orders (120 and 150). Despite the prominent spatial similarity of RSNs across model orders, some of the RSNs with significant differences at model order 100 and 120 do not show any differences in the equivalent RSNs at model order 150 and vice versa (see Figures 3, 4).

Table 3

SAD functional connectivity differences in the motor cortex at all estimated ICA model orders.

Model order	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
40	Broca's area BA44 R	0.5	21	18	76	50	18.14	76.62	50.33
	Broca's area BA45 R	0.6	25	18	76	50	18.44	77.9	47.59
60	Broca's area BA44 R	1.1	46	24	65	56	25.11	66.54	57.78
	Broca's area BA44 L	0.1	4	66	76	34	66.74	75.75	34
	Broca's area BA45 L	0.1	6	66	76	34	66.33	75.84	33.84
	Premotor cortex BA6 R	0.1	15	38	72	56	31.59	68.92	57.26
70	Broca's area BA44 R	1.9	80	26	66	56	20.54	72.62	53.47
	Broca's area BA45 R	1.4	61	19	76	54	18.32	75.62	51.87
	Premotor cortex BA6 R	0.5	62	38	72	56	31.63	66.91	55.7
80	Broca's area BA44 R	0.5	19	18	76	50	19.26	75.89	50.26
	Broca's area BA45 R	0.8	35	18	76	50	19.43	77.82	49.88
100	Broca's area BA44 R	1.5	64	24	65	56	22.16	69.29	52.32
	Broca's area BA45 R	0.9	38	22	70	50	20.21	73.25	49.24
	Premotor cortex BA6 R	0.8	109	34	72	60	32.41	69	61.64
	1ry motor cortex BA4a R	0.2	16	24	58	58	24.25	58.62	58.13
	1ry motor cortex BA4p R	0.3	14	24	58	58	25.78	56.79	57.15
	1ry somatosensory cortex BA3a R	0.3	10	24	58	58	26.5	56.1	56.6
	1ry somatosensory cortex BA3b R	0.2	13	24	58	58	24.38	57.85	58
120	Broca's area BA44 R	2.4	99	22	66	50	23.38	65.63	52
	Premotor cortex BA6 R	3.3	437	40	68	60	31.77	61.39	58.17
	1ry motor cortex BA4a R	1.7	116	26	58	54	24	58.06	56.72
	1ry motor cortex BA4p R	3.5	181	24	60	52	25.66	56.99	55.24
	1ry somatosensory cortex BA3a R	5.1	181	24	60	52	26.29	56.85	55.03
	1ry somatosensory cortex BA3b R	1.8	121	24	60	52	24.56	56.64	56.24
	1ry somatosensory cortex BA1 R	0.3	17	20	54	56	20.7	55.47	56.82
	1ry somatosensory cortex BA2 R	0.4	18	20	54	56	21.53	54.4	56.99
150	Broca's area BA44 R	0.1	3	22	68	50	22	68.66	51.66

The table describes one of the five motor resting-state networks detected at model order 40 and its subsequent (60–150) sub-networks. SAD significant connectivity differences (see Figure .

Table 7

SAD functional connectivity differences in the motor cortex at all estimated ICA model orders.

Model order	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
40	BA44 R	0.6	26	12	64	46	12.73	63.42	46.38
	1ry somatosensory cortex BA1 R	0.7	39	12	64	46	12.56	62.97	46.61
	1ry somatosensory cortex BA3b R	0.6	38	12	64	46	12.58	62.92	46.58
	1ry motor cortex BA4p R	0.6	29	12	64	46	12.83	62.83	46.82
	2ry somatosensory cortex OP4 R	1.1	35	12	64	46	12.54	63	46.54
60	BA44 R	0.5	21	12	64	46	13.05	63.47	46.85
	1ry somatosensory cortex BA1 R	0.8	46	12	64	46	13.56	61.82	49.77
	1ry somatosensory cortex BA3b R	0.6	43	12	64	46	13.62	61.69	49.67
	1ry motor cortex BA4a R	0.5	33	12	64	47	13.9	61.47	51.09
	1ry motor cortex BA4p R	0.7	39	12	64	46	13.58	61.76	49.23
	2ry somatosensory cortex OP4 R	1.0	30	12	64	46	12.9	62.97	47.23
70	BA44 R	1.0	41	12	62	46	12.64	64.09	46.51
	1ry somatosensory cortex BA1 R	2.4	135	12	62	46	13.49	61.25	49.88
	1ry somatosensory cortex BA2 R	1.1	56	12	62	46	14.27	59.43	50.66
	1ry somatosensory cortex BA3a R	2.2	76	14	64	46	15.31	59.97	49.15
	1ry somatosensory cortex BA3b R	2.2	144	12	62	46	14.43	60.66	50.08
	1ry motor cortex BA4a R	1.3	91	13	64	46	14	61.09	51.66
	1ry motor cortex BA4p R	2.3	120	12	62	46	14.87	60.38	50.49
	Premotor cortex BA6 R	0.8	99	12	63	46	13.82	61.63	51.02
	2ry somatosensory cortex OP4 R	2.5	79	12	62	46	12.67	62.23	47.02
80	BA44 R	0.6	26	12	64	46	12.77	63.65	46.61
	1ry somatosensory cortex BA1 R	1.1	60	12	64	46	13.63	61.47	50.93
	1ry somatosensory cortex BA3a R	0.5	16	13	63	46	13.44	62.38	46.81
	1ry somatosensory cortex BA3b R	0.8	56	12	64	46	13.68	61.4	50.78
	1ry motor cortex BA4a R	0.6	43	13	64	46	14.05	61	52.81
	1ry motor cortex BA4p R	0.9	49	12	64	46	13.83	61.27	50.99
	Premotor cortex BA6 R	0.4	48	12	64	46	13.83	61.54	51.85
	2ry somatosensory cortex OP4 R	1.1	34	12	64	46	12.71	63.24	46.76
100	1ry somatosensory cortex BA1 L	0.5	32	72	54	56	72.91	54.56	56.6
	1ry somatosensory cortex BA2 L	0.6	32	72	54	56	72.91	54.56	56.6
	1ry somatosensory cortex BA3b L	0.4	27	72	54	56	72.71	54.81	56.34
	1ry motor cortex BA4a L	0.3	25	72	54	56	72.72	54.88	56.56
	1ry motor cortex BA4p L	0.3	20	72	54	56	72.35	54.9	56.05
120	–	–	–	–	–	–	–	–	–
150	1ry somatosensory cortex BA1 L	0.6	38	72	54	56	73.16	56.9	54.12
	1ry somatosensory cortex BA2 L	0.6	36	72	54	56	72.86	56.38	54.4
	1ry somatosensory cortex BA3a L	0.9	36	72	60	48	72.03	59.42	49.3
	1ry somatosensory cortex BA3b L	0.8	59	72	54	56	72.4	57.75	52.06
	1ry motor cortex BA4a L	0.7	51	72	54	56	72.67	57.76	52.38
	1ry motor cortex BA4p L	0.9	51	72	54	56	72.16	58.17	51.29
	2ry somatosensory cortex OP4 L	0.6	17	72	58	50	71.12	58.06	49.88

The table describes one of the five motor resting-state networks detected at model order 40 and its subsequent (60–150) sub-networks. SAD significant connectivity differences (see Figure .

Discussion

Independent component analysis model order significantly influences the detected between-group differences in functional connectivity measured using the ICA dual regression approach. The total volume of significant between-group differences reaches maximum at model order 70, then further increases in ICA model order convey fewer between-group differences (Figure 1B). Also, the total volume of RSNs increased up to model order 70, after which it plateaus. Notably, the results revealed only increased functional connectivity in SAD at all estimated model orders facilitating inferences on the effects of model order selection.

The total volume of between-group differences and the total volume of RSNs showed a general trend of increase (Figure 1), suggesting a relationship between the spatial coverage of RSNs and that of between-group differences. However, there is a marked elevation between model orders 60 and 70 in the total volume of between-group differences without a corresponding increase in the total volume of RSNs. Also, while the total volume of RSNs plateaus at high model orders (70–150), the total volume of between-group differences gradually decreases. Moreover, at model orders 70 and 100, the fraction of significant voxels reaches maximum (see Figure S1 in Supplementary Material) which might be linked to the detected local maxima in the total volume of between-group differences. These findings might indicate possible contributions of other factors in the estimation of between-group differences.

Independent component analysis utilizes the entire spatial extent to estimate both large-scale RSNs and fine-grained RSNs by decomposing the functional data into components (ICA spatial modes) according to the estimated/selected model order (McKeown et al., 1998; Calhoun et al., 2001; Kiviniemi et al., 2003). Functional brain segmentation at model orders 20, 40, 60, 70, 80, 100, 120, and 150 yielded 27, 40, 47, 54, 70, 81, and 95 RSNs, respectively (Table 2). These RSN segmentations showed significant differences across model orders. Notably, at model order 20, 11 large-scale RSNs cover most of the brain cortex with relatively low z-score threshold of three. However, the most coherent core areas within these 11 RSNs at z-score > 5 do not involve all subcortical regions, parts of parietal and temporal cortices, or the cerebellum, c.f. Figure 2. On the other hand, fine-grained RSNs do cover the entire brain cortex as well as the cerebellum at the same threshold (z-score > 5).

The z-scores of the IC maps are non-linearly affected by model order as the calculation of z-score involves dividing by the standard deviation of the residual variance which decreases as a function of ICA model order. The RSN brain coverage perhaps could be adjusted accordingly by modifying the z-score threshold such that the coverage would be identical at all model orders. However, in the group-PICA dual regression approach applied here, the RSN brain coverage does not play a significant role in the results. The dual regression utilizes un-thresholded group-PICA maps in order to avoid double dipping of the data. Also, un-thresholded RSN volumes showed only 2.3 % difference in brain coverage (data not shown here). Moreover, the TFCE method (used to retain only significant between-group differences in dual regression) does not use such fixed thresholds either.

Although there is a high spatial consistency between RSNs at high model orders e.g., 70–150 (see Figures 3, 4), the detected between-group differences do not show the same degree of consistency (see Tables 3–7). It seems that some between-group differences are relatively more consistent across model orders than others. Each model order is a unique solution to the decomposition of the spatial brain data and subsequent model orders usually yield RSNs that do not share exactly the same brain regions. In spatial domain ICA, the mixing vector representing function in time is the one being modified when maximizing non-Gaussian distribution of spatial ICA modes. The time-domain mixing vector is likely to capture somewhat different temporal aspects of the data within each model order, particularly at high model orders that retain more of the original data for the actual ICA step. Although the core of RSNs remains very similar across model orders, between-group differences are uniquely estimated at each model order.

Our results showed that the total volume of RSN maps captured at high model orders exceeds that of large-scale RSNs at the same threshold p < 0.5 level by 13% (Figure 1, left). The exceeding volume could be explained by: (1) Filling in of the uncovered brain regions at low model orders, e.g., the cerebellum. (2) Partial overlaps within non-gray matter brain regions, e.g., white matter at a given threshold. In our analysis the PCA data reduction step prior to ICA was matched to the number of ICA model order. Thus, the non-gray matter brain regions shared by some RSNs are likely due to ICA capturing more of the data variance and/or detecting different T2*-weighted image effects and assigning them to different RSNs. In addition, at the high model orders (e.g., >100), it is sometimes difficult to differentiate RSN and non-RSN sources due to the same reason of increased variance capture.

Multilevel exploration of functional brain connectivity

Measuring functional brain connectivity at different hierarchical levels using ICA might provide an advantage in modeling the systematic effects of pathology in brain disorders. Probably, diagnosing neuropsychiatric disorders might require selecting the optimal ICA model order according to the function under investigation, not only the number of RSNs or the volume of between-group differences. Also, information regarding the underlying functionality of fine-grained sub-networks shall play an essential role in determining the hierarchical level needed to profile each disorder. Additionally, in some brain disorders functional connectivity abnormalities might be detected among both large-scale and fine-grained RSNs, while in other disorders abnormalities may lay either among large-scale or fine-grained RSNs. Therefore, there may be no definite model order for detecting aberrant functional connectivity for all diseases.

As we mentioned earlier, SAD was used only as an example of brain disease in this study. The main focus of this study was to explore if disease sensitivity could be optimized by exploring functional connectivity differences at multiple ICA model orders. Therefore, detailed interpretation of SAD-related connectivity differences is beyond the scope of this paper. More interestingly however, our results showed SAD-related increased functional connectivity in RSNs including brain regions that have been reported to involve increased metabolism in SAD, e.g., in the frontal, orbitofrontal, right parietal, and left middle temporal cortices, and in the right caudate: whereas brain regions of reduced metabolism were also recorded, e.g., the right middle temporal cortex and left parietal cortex (Cohen et al., 1992). Increased functional connectivity might be associated with increased underlying cerebral metabolism. In addition, other brain areas not yet associated with SAD were also detected, e.g., the PCC, precuneus, and the visual cortex. Future work aims at detailed interpretation and discussion of SAD-related functional connectivity differences detected in this study on neuropsychiatric and neurofunctional bases.

Modularity demonstrated in a wide range of complex systems enables the identification of the optimal partition of a network. There is strong evidence for the modular organization of human brain networks (Chen et al., 2008; Bullmore and Sporns 2009; Ferrarini et al., 2009; Meunier et al., 2009a, 2009b; Valencia et al., 2009). Salvador et al. (2005) proposed partial correlation methods to build group-level RSNs based upon an averaged similarity matrix across subjects, while van den Heuvel et al. (2008) have implemented an individual clustering adjacency matrix. Also, Ferrarini et al. (2009) have described the hierarchical modular structure of RSNs using a partial correlation analysis. Meunier et al. (2009b) functionally segmented the brain into eight large modules, each compromising more than 10 nodes at the highest level of the hierarchy, and 57 sub-modules at the lowest level. Recently, Bellec et al. (2010) detected a minimum of seven RSNs using bootstrap analysis of stable clusters and suggested that there might be many numbers of clusters (where local maxima of clustering stability could be spotted). Our ICA segmentations are consistent with these previous results that highlight the modular segmentation of the human brain. Nevertheless, a more detailed investigation driven by similarity between the ICA functional hierarchy of RSNs and the above findings using different methodologies for assessing resting-state brain activity would be beneficial.

The potential utility of group based ICA approaches has been demonstrated by the increasing number of studies examining clinical populations. Large-scale networks (i.e., low model order components) are compact and easy to identify, including networks such as the visual, auditory, sensorimotor, etc. Notably, ICA shows tendency toward “splitting” of large-scale networks into sub-networks by increasing the number of estimated ICs. Previously, we showed that the use of high model order ICA is an important and useful tool for obtaining detailed parcellation of the functional hierarchy of neural sub-systems (Abou Elseoud et al., 2010). It is important to notice that different ICA model orders describe different functional connectivity and hierarchical aspects. Our results showed that a large-scale network with significant differences in connectivity splits into sub-networks. Some of these sub-networks still show connectivity differences, while the rest show normal connectivity, i.e., some of the sensorimotor and visual sub-networks show no changes in connectivity at higher model orders. This finding indicates that although the branched sub-networks descended from a common large-scale network, they are independent from each other at the fine-grained levels of the hierarchy.

It is also important to consider the neurofunctional meaning of such fine-grained sub-networks when selecting ICA model order. At present, assigning every sub-network to a neuronal sub-functionality appears to be difficult. A number of studies showed conclusions linking networks’ splitting and sub-functionalities. Seeley et al. (2007) identified such a split in a network suggested to be involved in executive function, revealing separate purported “control” and “salience” networks. Earlier, Seifritz et al. (2002) showed ICA of the auditory cortex that yielded a unique spatial map and a temporal pattern, suggesting the presence of two concurrent temporally independent processes. Then by temporally decomposing the signal into temporal ICs, transient and sustained components of the human auditory cortex were separated. Additionally, high model orders can be used to detect individual subcortical components belonging to the thalamus and the basal ganglia, which are not detectable using low model order (Ystad et al., 2010). Moreover, on the millimeter-scale, marked differences in the patterns of functional connectivity can be detected (Margulies et al., 2007). Future studies are needed in order to link functional segmentation to anatomical parcellation and to fill in the gap between fine-grained networks and sub-functionality.

At high model orders there is a large increase in statistical tests performed compared to low model orders due to the increased number of functional RSN sources. The presently applied TFCE correction method corrects the results at brain voxel level, but does not adjust for the risk (false positives) induced by increasing the number of RSNs. Currently we are planning statistical methods that could correct for the type I error also at the level of RSN in combination with dual regression. Therefore, we performed corrections only at individual IC level. In general, our results at high model orders, particularly model order 70 are in line with results at low model orders. Moreover, the total volume of significant between-group differences showed no significant increase after model order 70. However, it would be beneficial to develop a methodology for correcting for multiple comparisons simultaneously within and between RSN maps.

Functional hierarchy and disease sensitivity

Neuropsychiatric disorders (e.g., depression, dementia, schizophrenia, Alzheimer's, etc.) might occur at either the micro-level or macro-level, particularly at the initial stages, and may extend to both levels later on (Stoffers et al., 2008; Krishnan et al., 2009; Zhang et al., 2010). Some disorders target specific brain systems or regions and spare others within the same anatomical structure, consequently affecting some functions and preserving others. Importantly, the choice of ICA model order when the ultimate goal is detecting between-group difference should be considered cautiously. Our results showed a small number of brain voxels with significant between-group differences at low model orders (Figure 1, right). Moreover, large-scale RSNs at lower model orders suffer from the joining of multiple areas known to be functionally independent (Abou Elseoud et al. 2010; see Figure 3 and Table 3, at model order 20, one RSN included both motor and auditory cortices). Also, large-scale RSNs are unable to capture subcortical brain regions as separate RSNs (Abou Elseoud et al. 2010; Ystad et al., 2010; see Figure 2). On the contrary, at high model orders the total volume of between-group differences reaches maximum. However, there is an increased risk of false positives due to an increased number of components at higher model orders. Therefore, as discussed earlier it might be that applying different statistical methods for correcting results at high model orders could provide a reasonable solution.

We emphasize that results presented here apply to our 1.5 Tesla group-ICA setting with the present imaging parameters and data pre-processing. Additionally, different optimal model orders might be found when higher field strengths and higher resolutions are used. However, as clinical applications require analysis of individual fMRI datasets, an additional challenge remains, concerning the matter of how one can detect subject-specific RSNs robustly at fine-grained levels. High field magnets, e.g., 3 or 7 Tesla, might offer an alternative which could provide individual subject datasets with adequate quality needed for detecting subject-specific RSNs.

Conclusion

In this paper we show that detected disease-related differences in functional connectivity alter as a function of ICA model order. The results showed only increased functional connectivity at all estimated model orders. Our findings suggest that multilevel ICA exploration of functional connectivity enables optimization of sensitivity to brain disorders. The volume of between-group differences altered significantly as a function of ICA model order reaching maximum at model order 70 (which seems to be an optimal level that conveys the largest between-group difference). Our results show that fine-grained RSNs enable better detection of detailed disease-related functional connectivity changes. However, high model orders show an increased risk of false positives that needs to be overcome.

Conflict of Interest Statement

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

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/systems_neuroscience/10.3389/fnsys.2011.00037/abstract

Table 4

SAD functional connectivity differences in the motor cortex at all estimated ICA model orders.

Model order	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
40	–	0	0	–	–	–	–	–	–
60	–	0	0	–	–	–	–	–	–
70	1ry motor cortex BA4a L	0.3	24	62	56	54	63.3	55.63	55.16
	Premotor cortex BA6 R	0.1	8	34	56	60	40	62.49	59
	Premotor cortex BA6 L	1.1	137	50	66	54	49.45	67.51	58.75
	1ry motor cortex BA4p L	0.6	33	62	56	54	63.9	55.46	54.72
	1ry somatosensory cortex BA3a L	0.8	33	62	56	54	63.9	55.46	54.72
	1ry somatosensory cortex BA3b L	0.4	32	62	56	54	63.96	55.38	54.75
	Anterior intra-parietal sulcus hIP1 L	1.1	38	62	32	50	62.86	31.86	50.95
	Anterior intra-parietal sulcus hIP2 L	0.7	20	62	32	50	62.54	32.59	51.15
80	1ry motor cortex BA4a R	0.1	6	45	56	60	44.5	55.5	60
	1ry motor cortex BA4a L	0.4	29	46	56	60	46.62	55.59	59.62
	Premotor cortex BA6 R	0.3	38	45	68	58	43.84	65.51	58.29
	Premotor cortex BA6 L	0.9	107	46	68	58	47.11	64.1	59.14
100	1ry motor cortex BA4a L	7.1	546	66	53	56	67.46	57.44	56.51
	Premotor cortex BA6 L	3.0	369	66	55	58	67.95	58.34	57.54
	1ry somatosensory cortex BA1 L	3.9	258	66	51	56	69.13	56.68	57.79
	1ry motor cortex BA4p L	11.4	674	62	52	54	66.41	55.92	55.59
	1ry somatosensory cortex BA2 L	5.1	294	65	52	54	68.33	54.55	56.37
	1ry somatosensory cortex BA3a L	12.9	534	62	52	54	65.31	55.22	54.56
	1ry somatosensory cortex BA3b L	8.6	660	62	52	54	66.76	55.67	55.91
	Anterior intra-parietal sulcus hIP2 L	1.9	56	66	52	55	66.5	49.97	52
	2ry somatosensory cortex OP4 L	0.9	26	69	56	52	70.42	56.69	56.08
120	1ry motor cortex BA4a L	1.8	138	62	56	54	64.48	55.93	55.24
	Premotor cortex BA6 L	0.5	56	60	56	56	62.71	56.1	56.67
	1ry motor cortex BA4p L	4.9	292	64	52	54	63.91	54.86	54.72
	1ry somatosensory cortex BA2 L	1.7	98	65	52	54	66.29	52.29	54.22
	1ry somatosensory cortex BA3a L	8.1	336	64	52	54	63.85	54.49	54.46
	1ry somatosensory cortex BA3b L	3.6	274	64	52	54	64.37	54.29	54.64
	Anterior intra-parietal sulcus hIP2 L	1.2	35	64	50	53	65.34	49.37	53.85
	2ry somatosensory cortex OP2 L	2.4	20	64	48	49	62.55	47.84	48.3
	2ry somatosensory cortex OP1 L	1.2	43	64	48	50	63.46	48.64	49.25
	1ry auditory cortex TE1.1 L	1.5	22	64	48	49	62.18	47.59	48.27
150	1ry motor cortex BA4a L	2.0	156	62	56	54	64.97	55.02	56.18
	Premotor cortex BA6 L	0.5	57	60	54	56	63.9	55.82	57.53
	1ry motor cortex BA4p L	4.6	272	64	52	53	64.19	54.28	55.31
	1ry somatosensory cortex BA2 L	1.9	109	65	52	54	65.96	52.78	55.57
	1ry somatosensory cortex BA3a L	7.2	298	64	51	52	64.04	53.93	54.92
	1ry somatosensory cortex BA3b L	3.6	274	64	52	52	64.62	54.05	55.39
	Anterior intra-parietal sulcus hIP2 L	0.6	18	64	50	53	64.67	49.11	52.99

The table describes one of the five motor resting-state networks detected at model order 40 and its subsequent (60–150) sub-networks. SAD significant connectivity differences (see Figure .

Table 5

SAD functional connectivity differences in the motor cortex at all estimated ICA model orders.

Model order	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
40	Anterior intra-parietal sulcus hIP2 R	0.6	18	30	50	58	28.9	49.83	58.77
	1ry motor cortex BA4a R	4.6	320	30	50	64	36.76	46.84	62.28
	1ry motor cortex BA4a L	2.7	203	62	56	60	61.16	52.18	60.19
	1ry motor cortex BA4p R	6.4	335	30	52	56	32.15	49.55	60.91
	1ry motor cortex BA4p L	4.3	252	62	56	60	60.85	51.41	59.74
	1ry somatosensory cortex BA1 R	0.9	51	30	50	64	30.65	49.33	65.49
	1ry somatosensory cortex BA1 L	0.9	63	66	54	58	64.86	51.55	59.38
	1ry somatosensory cortex BA2 R	0.6	29	29	50	60	29.83	48.62	59.38
	1ry somatosensory cortex BA2 L	2.1	123	65	54	58	61.33	48.24	59.06
	1ry somatosensory cortex BA3a R	9.1	321	30	52	56	32.12	49.65	60.11
	1ry somatosensory cortex BA3a L	5.3	219	62	56	59	59.06	49.21	59.27
	1ry somatosensory cortex BA3b R	4.1	275	30	52	57	30.83	50.17	60.76
	1ry somatosensory cortex BA3b L	3.1	238	63	56	60	61.53	51.48	59.37
	Premotor cortex BA6 R	1.3	176	30	52	56	31.78	50.89	60.98
	Premotor cortex BA6 L	1.1	132	62	56	60	62.46	53.98	60.25
60	–	–	–	–	–	–	–	–	–
70	Anterior intra-parietal sulcus hIP1 R	1.8	68	24	38	60	22.76	39.57	59.81
	Anterior intra-parietal sulcus hIP1 L	2.8	94	60	36	58	60.86	37	59.23
	Anterior intra-parietal sulcus hIP2 R	2.4	68	24	38	60	23.23	40.59	59.96
	Anterior intra-parietal sulcus hIP2 L	1.4	41	60	38	58	62.5	39.09	59.67
	1ry motor cortex BA4a R	0.8	56	38	53	60	38.95	50.37	61.87
	1ry motor cortex BA4a L	1.0	73	60	50	61	59.57	50.04	61.32
	1ry motor cortex BA4p R	1.0	53	36	46	61	35.57	48.16	61.45
	1ry motor cortex BA4p L	2.5	148	60	50	60	59.36	48.55	60.68
	1ry somatosensory cortex BA1 L	0.6	39	60	50	61	60.84	46.75	61.03
	1ry somatosensory cortex BA2 R	0.6	28	24	42	60	27.58	45.1	60.71
	1ry somatosensory cortex BA2 L	2.2	125	60	48	62	60.57	44.53	60.31
	1ry somatosensory cortex BA3a R	2.0	69	36	46	60	35.87	47.54	61.14
	1ry somatosensory cortex BA3a L	3.7	154	60	50	60	59.36	48.34	60.61
	1ry somatosensory cortex BA3b L	1.8	138	60	50	60	59.43	48.58	60.79
	Premotor cortex BA6 R	0.6	77	38	46	60	38.8	50.64	61.56
80	–	–	–	–	–	–	–	–	–
100	–	–	–	–	–	–	–	–	–
120	–	–	–	–	–	–	–	–	–
150	–	–	–	–	–	–	–	–	–

The table describes one of the five motor resting-state networks detected at model order 40 and its subsequent (60–150) sub-networks. SAD significant connectivity differences (see Figure .

Table 6

SAD functional connectivity differences in the motor cortex at all estimated ICA model orders.

Model order	Anatomical region	Percent of overlap	Non-zero voxels	Max z-score			Center of mass
				X	Y	Z	X	Y	Z
60 A	Anterior intra-parietal sulcus hIP1 R	0.2	8	21	44	64	22.24	43.51	63.12
	Anterior intra-parietal sulcus hIP2 R	0.7	21	20	44	64	22.04	44.76	62.95
	1ry somatosensory cortex BA1 R	0.2	13	20	44	64	21.69	45.3	63.23
	1ry somatosensory cortex BA2 R	0.4	21	20	44	64	22.04	44.76	62.95
60 B	Anterior intra-parietal sulcus hIP1 L	0.9	32	62	42	58	64	40.73	61.08
	Anterior intra-parietal sulcus hIP2 L	2.6	75	68	48	60	66.01	43.55	61.41
	1ry motor cortex BA4a L	3.2	248	68	50	60	65.78	50.9	62.36
	1ry motor cortex BA4p L	4.3	255	70	52	60	65.38	50.17	61.97
	1ry somatosensory cortex BA1 L	5.3	354	68	48	60	66.19	48.73	62.18
	1ry somatosensory cortex BA2 L	6.3	363	68	48	60	66.05	47.8	61.69
	1ry somatosensory cortex BA3a L	3.2	134	68	49	60	64.9	48.41	61.26
	1ry somatosensory cortex BA3b L	4.6	351	68	48	60	65.91	49.49	61.99
	Premotor cortex BA6 L	0.7	82	62	54	62	64.36	53.82	63.12
60 C	–	–	–	–	–	–	–	–	–
70 A	Anterior intra-parietal sulcus hIP1 R	3.1	119	20	44	63	24.1	43.28	61.04
	Anterior intra-parietal sulcus hIP2 R	13.7	385	20	44	64	21.38	46.33	60.28
	1ry somatosensory cortex BA1 R	7.1	404	20	44	64	24.98	48.03	63.26
	1ry somatosensory cortex BA2 R	9.4	462	20	44	64	22.87	46.49	61.53
	1ry somatosensory cortex BA3a R	5.1	179	26	46	60	28.86	48.97	62.75
	1ry somatosensory cortex BA3b R	5.5	363	25	46	60	26.9	49.38	63.07
	2ry somatosensory cortex OP1 R	3.2	91	18	46	56	18.89	48.31	55.75
	Premotor cortex BA6 R	2.1	272	30	50	64	28.52	51.57	64.47
	1ry motor cortex BA4a R	4.3	301	30	50	64	28.2	51.07	64.41
	1ry motor cortex BA4p R	5.8	300	26	46	60	28.94	49.91	64.15
70 B	Anterior intra-parietal sulcus hIP1 L	1.4	47	65	43	60	64.78	39.72	61.03
	Anterior intra-parietal sulcus hIP2 L	3.5	102	66	44	60	65.93	42.4	61.24
	1ry motor cortex BA4a L	2.9	223	66	46	63	64.78	48.73	63.43
	1ry motor cortex BA4p L	4.2	246	64	44	62	64.16	47.93	62.76
	1ry somatosensory cortex BA1 L	5.2	346	66	44	60	65.46	47.31	62.92
	1ry somatosensory cortex BA2 L	6.9	399	66	44	60	65.15	46.5	62.25
	1ry somatosensory cortex BA3a L	4.0	166	66	46	60	64.14	47.5	61.52
	1ry somatosensory cortex BA3b L	4.5	345	66	46	60	64.9	47.72	62.73
	Premotor cortex BA6 L	0.4	47	64	48	66	63.34	51.22	65.79
70 C	–	–	–	–	–	–	–	–	–
80 A	Anterior intra-parietal sulcus hIP1 R	0.5	18	21	44	64	20.94	43.67	62.78
	Anterior intra-parietal sulcus hIP2 R	2.0	57	20	44	64	21.19	44.94	63.09
	1ry somatosensory cortex BA1 R	1.0	55	20	44	64	23.53	47.41	63.62
	1ry somatosensory cortex BA2 R	1.2	61	20	44	64	21.14	45.14	63.17
	1ry somatosensory cortex BA3a R	0.3	9	26	46	60	28.43	48.77	63.21
	1ry somatosensory cortex BA3b R	0.5	30	30	50	66	26.63	49.59	63.4
	1ry motor cortex BA4a R	0.3	20	30	50	66	27.91	50.99	64.55
	1ry motor cortex BA4p R	0.4	19	30	50	66	29.47	49.42	64.89
80 B	Anterior intra-parietal sulcus hIP1 L	0.9	29	62	44	58	64.96	39.52	60.96
	Anterior intra-parietal sulcus hIP2 L	2.6	74	66	44	60	66.04	43.1	61.12
	1ry motor cortex BA4a L	1.6	122	66	47	62	65.93	48.39	62.56
	1ry motor cortex BA4p L	2.2	128	66	46	62	65.23	47.73	62.13
	1ry somatosensory cortex BA1 L	3.7	248	66	44	60	66.06	46.63	62.35
	1ry somatosensory cortex BA2 L	4.9	283	66	44	60	65.77	46.33	62.03
	1ry somatosensory cortex BA3a L	2.3	97	66	46	60	65.16	47.74	61.29
	1ry somatosensory cortex BA3b L	2.9	221	66	46	60	65.79	47.13	62.33
80 C	–	–	–	–	–	–	–	–	–
100 A	Anterior intra-parietal sulcus hIP1 R	0.4	14	20	44	63	20.93	43.79	62.58
	Anterior intra-parietal sulcus hIP2 R	0.8	23	20	44	64	20.69	44.18	63.09
	1ry somatosensory cortex BA1 R	0.2	10	20	44	64	20.6	44.79	63.9
	1ry somatosensory cortex BA2 R	0.5	23	20	44	64	20.69	44.18	63.09
100 B	Anterior intra-parietal sulcus hIP1 L	4.5	154	64	41	62	64.71	38.26	61.63
	Anterior intra-parietal sulcus hIP2 L	7.0	202	66	44	60	65	40.74	61.64
	1ry motor cortex BA4a L	2.6	197	66	47	62	65.12	49.15	63.1
	1ry motor cortex BA4p L	3.3	193	66	46	62	64.83	48.66	62.55
	1ry somatosensory cortex BA1 L	6.2	414	66	44	60	65.02	45.71	63.21
	1ry somatosensory cortex BA2 L	8.3	483	66	44	60	64.88	44.89	62.61
	1ry somatosensory cortex BA3a L	3.0	123	66	47	62	64.96	48.29	61.37
	1ry somatosensory cortex BA3b L	4.0	303	66	45	62	65.2	47.64	62.79
	Premotor cortex BA6 L	0.5	58	64	52	62	62.86	54.07	65.91
100 C	Anterior intra-parietal sulcus hIP1 L	0.9	31	62	44	58	65.25	42.91	58.83
	Anterior intra-parietal sulcus hIP2 L	1.2	36	64	44	58	66.19	43.14	59.3
	1ry motor cortex BA4p L	0.7	40	70	56	48	70.06	54.98	49.44
	1ry somatosensory cortex BA1 L	0.3	18	70	53	50	68.22	47.4	57.03
	1ry somatosensory cortex BA2 L	1.2	72	70	56	49	67.58	47.98	55.46
	1ry somatosensory cortex BA3a L	1.5	64	70	56	48	69.68	54.18	49.5
	1ry somatosensory cortex BA3b L	1.0	75	70	56	48	69.85	53.68	49.73
	2ry somatosensory cortex OP1 L	2.1	76	70	56	48	70.07	53.7	49.15
	2ry somatosensory cortex OP2 L	1.7	14	69	54	48	68.43	54.86	48.57
	2ry somatosensory cortex OP3 L	1.1	18	70	56	48	69.01	54.89	48.28
	2ry somatosensory cortex OP4 L	2.7	81	70	56	48	70.04	53.94	49.1
120 A	Anterior intra-parietal sulcus hIP1 R	0.3	10	20	48	56	20.4	44.61	61.29
	Anterior intra-parietal sulcus hIP2 R	1.5	41	18	48	56	19.36	47.23	58.28
	1ry somatosensory cortex BA1 R	0.8	47	18	50	56	23.48	48.84	61.74
	1ry somatosensory cortex BA2 R	1.2	61	18	49	56	20.64	48.73	58.38
	1ry somatosensory cortex BA3a R	2.4	84	28	54	60	27.5	53.34	58.86
	1ry somatosensory cortex BA3b R	1.3	89	28	54	60	27.19	52.4	59.78
	2ry somatosensory cortex OP1 R	1.8	50	18	48	56	18.5	49.86	55.46
	Premotor cortex BA6 R	0.5	72	28	54	60	28.36	52.18	61.99
	1ry motor cortex BA4a R	1.1	74	28	54	60	28.2	52.21	61.96
	1ry motor cortex BA4p R	1.9	99	28	54	60	27.84	52.81	60.04
120 B	Anterior intra-parietal sulcus hIP1 L	0.9	29	62	40	60	62.79	40.11	60.45
	Anterior intra-parietal sulcus hIP2 L	1.2	36	62	40	60	62.94	40.5	60.58
	1ry motor cortex BA4a L	0.1	5	56	46	60	56	47.2	59.6
	1ry motor cortex BA4p L	0.5	31	56	44	60	57.2	44.86	60.2
	1ry somatosensory cortex BA1 L	0.9	59	61	40	60	59.35	41.48	61.12
	1ry somatosensory cortex BA2 L	1.8	106	62	40	60	59.8	41.77	60.75
	1ry somatosensory cortex BA3a L	1.4	56	58	42	60	56.97	44.13	60.11
	1ry somatosensory cortex BA3b L	0.7	52	56	42	60	57.2	43.36	60.73
120 C	Anterior intra-parietal sulcus hIP1 L	0.2	7	60	44	58	63.41	43.72	58.42
	Anterior intra-parietal sulcus hIP2 L	1.4	42	65	43	62	66.44	44.94	59.02
	1ry motor cortex BA4a L	0.2	15	64	44	64	64.67	46.26	62.6
	1ry motor cortex BA4p L	1.0	57	64	44	62	63	45.72	60.77
	1ry somatosensory cortex BA1 L	1.1	71	64	44	62	65.36	44.71	61.89
	1ry somatosensory cortex BA2 L	2.2	125	64	44	62	64.49	45.1	60.38
	1ry somatosensory cortex BA3a L	1.3	52	64	44	62	62.58	45.67	60.42
	1ry somatosensory cortex BA3b L	1.1	82	64	44	62	63.72	45.38	61.07
150 A	Anterior intra-parietal sulcus hIP1 R	0.2	9	20	48	56	22.86	48	56.33
	Anterior intra-parietal sulcus hIP2 R	3.7	103	18	49	56	22.97	48.9	57.04
	1ry somatosensory cortex BA1 R	2.4	138	18	50	56	25.79	49.55	63.01
	1ry somatosensory cortex BA2 R	3.2	156	18	50	56	22.63	49.44	58.04
	1ry somatosensory cortex BA3a R	4.7	167	30	50	64	27	51.28	58.67
	1ry somatosensory cortex BA3b R	3.7	246	30	50	64	26.62	50.79	60.75
	2ry somatosensory cortex OP1 R	2.4	69	18	50	56	20.05	50.41	54.6
	Premotor cortex BA6 R	1.2	159	30	50	64	29.26	51.25	64.68
	1ry motor cortex BA4a R	2.5	173	30	50	64	28.96	51.01	64.85
	1ry motor cortex BA4p R	4.3	225	30	50	64	28.32	51.24	62.47
	2ry somatosensory cortex OP2 R	0.7	9	26	56	50	25.33	53.44	51
	2ry somatosensory cortex OP3 R	1.0	22	24	58	50	25.22	56.73	50.49
150 B	–	–	–	–	–	–	–	–	–
150 C	Anterior intra-parietal sulcus hIP1 L	0.1	3	62	46	54	62	45	54
	Anterior intra-parietal sulcus hIP2 L	0.3	8	67	45	62	65.53	44.88	59.05
	1ry motor cortex BA4a L	0.4	33	66	47	62	67.05	48.84	62.51
	1ry motor cortex BA4p L	0.6	36	64	44	62	66.18	47.79	62.46
	1ry somatosensory cortex BA1 L	1.0	64	64	44	62	66.6	47.1	62.48
	1ry somatosensory cortex BA2 L	1.2	67	64	44	62	66.48	46.99	62.41
	1ry somatosensory cortex BA3a L	0.6	26	64	44	62	65.69	47.46	62.3
	1ry somatosensory cortex BA3b L	0.8	60	64	44	62	66.49	47.3	62.51

The table describes one of the five motor resting-state networks detected at model order 40 and its subsequent (60–150) sub-networks. SAD significant connectivity differences (see Figure .