Multimodel Order Independent Component Analysis: A Data-Driven Method for Evaluating Brain Functional Network Connectivity Within and Between Multiple Spatial Scales

Abstract

Background:

While functional connectivity is widely studied, there has been little work studying functional connectivity at different spatial scales. Likewise, the relationship of functional connectivity between spatial scales is unknown.

Methods:

We proposed an independent component analysis (ICA)-based approach to capture information at multiple-model orders (component numbers), and to evaluate functional network connectivity (FNC) both within and between model orders. We evaluated the approach by studying group differences in the context of a study of resting-state functional magnetic resonance imaging (rsfMRI) data collected from schizophrenia (SZ) individuals and healthy controls (HC). The predictive ability of FNC at multiple spatial scales was assessed using support vector machine-based classification.

Results:

In addition to consistent predictive patterns at both multiple-model orders and single-model orders, unique predictive information was seen at multiple-model orders and in the interaction between model orders. We observed that the FNC between model orders 25 and 50 maintained the highest predictive information between HC and SZ. Results highlighted the predictive ability of the somatomotor and visual domains both within and between model orders compared with other functional domains. Also, subcortical-somatomotor, temporal-somatomotor, and temporal-subcortical FNCs had relatively high weights in predicting SZ.

Conclusions:

In sum, multimodel order ICA provides a more comprehensive way to study FNC, produces meaningful and interesting results, which are applicable to future studies. We shared the spatial templates from this work at different model orders to provide a reference for the community, which can be leveraged in regression-based or fully automated (spatially constrained) ICA approaches.

Impact statement

Multimodel order independent component analysis (ICA) provides a comprehensive way to study brain functional network connectivity within and between multiple spatial scales, highlighting findings that would have been ignored in single-model order analysis. This work expands upon and adds to the relatively new literature on resting functional magnetic resonance imaging-based classification and prediction. Results highlighted the differentiating power of specific intrinsic connectivity networks on classifying brain disorders of schizophrenia patients and healthy participants, at different spatial scales. The spatial templates from this work provide a reference for the community, which can be leveraged in regression-based or fully automated ICA approaches.

Introduction

Brain activity reveals exquisite coordination between spatial scales, from local microcircuits to brain-wide networks. Brain activity measured by noninvasive functional brain imaging techniques is typically assumed to be generated on the cortical surface. The spatial extent of activity on the cortex obtained experimentally from neuroimaging modalities and models varies widely (Perdue and Diamond, 2013). Understanding the brain requires an integrated understanding of the different scales of the spatial organization of the brain. Studies that utilize measurements between spatial scales promise to increase our understanding of brain function by tracking sensory, motor, and cognitive variables as they evolve through local microcircuits and across brain-wide networks (Lewis et al., 2015).

Whole-brain functional connectivity can be studied by calculating the functional network connectivity (FNC) (Allen et al., 2011) between intrinsic connectivity networks (ICNs); that is, the functional magnetic resonance imaging (fMRI) independent components (IC) resulted from group independent component analysis (ICA). To date, most previous research has focused on FNC within one specific model order (i.e., number of components), ignoring the importance of capturing functional information at different levels of spatial granularity as well as the between-order information. A wide number of studies have focused on FNC in brain health and various disorders, all at a single-model order. In particular, multiple studies have highlighted significant FNC differences in studies of schizophrenia (SZ), for example, increased FNC in SZ among frontal and temporal networks (Calhoun et al., 2009), within the default mode networks (Salvador et al., 2010; Whitfield-Gabrieli et al., 2009), and decreased FNC between temporal and parietal networks (Jafri et al., 2008). While it is well known that there are functional changes associated with SZ, it is unclear to what degree these are linked to the choice of a specific model order (IC numbers) when performing FNC analysis.

Considering the brain as a functional hierarchy (Iraji et al., 2019), it is clear that functional interactions between functional sources can occur between different spatial scales, and the lack of studies that evaluate functional connectivity between spatial scales represents a gap in the field. Here, we propose to study FNC within and between four different functional hierarchy levels (model order = 25, 50, 75, and 100), thus providing new insight to understand brain functional connectivity. In this study, our goal is to compare the important information obtained at different model orders by evaluating the multiorder FNCs using healthy control (HC)–SZ group comparison, as well as testing the classification power of FNC in discriminating SZs based on machine learning-based approaches (Anderson and Cohen, 2013; Pariyadath et al., 2014). We seek to identify the FNC-based features that are predictive of SZ from controls through support vector machine (SVM)-based classifications. Specifically, we evaluated if the HC–SZ group differences primarily occur at one specific model order or the between-model orders, and which order yields the highest classification power between HC and SZ. The ICNs template at four model orders was also made available to the community for their future use.

Methods

Dataset and preprocessing

Our dataset was combined across three separate studies, one with seven sites (FBIRN: Functional Imaging Biomedical Informatics Research Network), one with three sites (MPRC: Maryland Psychiatric Research Center), and one single site (COBRE: Center for Biomedical Research Excellence). We extracted a subset of subjects from the three datasets satisfying the following criteria: (1) data of individuals with typical control or SZ diagnosis; (2) data with high-quality registration to echo-planar imaging (EPI) template; (3) head motion transition of <3° rotations and 3-mm translations in all directions. The mean framewise displacement (meanFD) among selected subjects was average ± standard deviation = 0.1778 ± 0.1228; min–max = 0.0355–0.9441. This resulted in 827 individuals, including 477 subjects (age: 38.76 ± 13.39, females: 213, males: 264) of HCs and 350 SZ individuals (age: 38.70 ± 13.14, females: 96, males: 254). The parameters of resting-state fMRI (rsfMRI) for the FBIRN data were the same across all sites, with a standard gradient EPI sequence, repetition time (TR)/echo time (TE) = 2000/30 ms, and a total of 162 volumes. For COBRE data, rsfMRI images were acquired using a standard EPI sequence with TR/TE = 2000/29 ms and 149 volumes. The MPRC datasets were acquired using a standard EPI sequence at three sites, including Siemens 3-Tesla Siemens Allegra scanner (TR/TE = 2000/27 ms, voxel spacing size = 3.44 × 3.44 × 4 mm, field of view [FOV] = 220 × 220 mm, and 150 volumes), 3-Tesla Siemens Trio scanner (TR/TE = 2210/30 ms, voxel spacing size = 3.44 × 3.44 × 4 mm, FOV = 220 × 220 mm, and 140 volumes), and 3-Tesla Siemens Tim Trio scanner (TR/TE = 2000/30 ms, voxel spacing size = 1.72 × 1.72 × 4 mm, FOV = 220 × 220 mm, and 444 volumes).

We were following the same rsfMRI data preprocessing procedures as reported in a recent study (Iraji et al., 2021). First of all, for the purposes of magnetization equilibrium, the first five volumes were discarded. We performed rigid motion correction to correct the head motion for each subject during scanning. We then applied slice-time correction to deal with temporal misalignment. And the next step, the rsfMRI data of each subject were registered to a Montreal Neurological Institute EPI template, resampled to 3 mm³ isotropic voxels, and smoothed using a Gaussian kernel with a 6 mm full width at half-maximum (6 mm). The voxel time courses were then z-scored (variance normalized). The preprocessing procedures were mainly performed by using the statistical parametric mapping (SPM12) toolbox.

Group ICA analysis

ICA provides adaptive, overlapping, networks at different scales; we and others have shown this has substantial advantages over fixed ROIs (Du et al., 2020; Yu et al., 2018). The main advantages of ICA analysis include its ability to adapt to individual subjects, and the ICA is optimized to include temporally coherent voxels. The overlapping networks captured by ICA analysis provide a natural decomposition of the data. In ICA analysis, the selection of model order (i.e., number of components to be extracted) can effectively define spatial scale of ICNs, and therefore has crucial effects on brain functional network analysis. This makes ICA a great tool to obtain subject-specific ICNs at different spatial scales. In another word, we can study brain function within and between different spatial scales by using ICA with different model orders. Low model order ICAs produce large-scale spatially distributed ICNs, such as the default mode network (Beckmann et al., 2005; Calhoun et al., 2008; Damoiseaux et al., 2008), whereas high model orders produce spatially granular ICNs, including multiple more focal ICNs instead of one large-scale default mode network (Allen et al., 2011). Given these advantages, we tend to study brain function at multiple scales using ICA.

The ICA analysis was performed using the GIFT software (Calhoun and Adali, 2012; Calhoun et al., 2001). In our study, ICA was performed at four model orders: 25, 50, 75, and 100. Before the ICA, subject-specific principal components analysis (PCA) was first performed on the datasets to normalize the data. The subject-level principal components were then concatenated together across the time dimensions, and the group-level PCA was applied to concatenate the subject-level principal components. The group-level principal components that explained the maximum variance were selected as the input of ICA to perform group-level ICA (Calhoun et al., 2001). We used the infomax algorithm and controlled for stochastic variability by using ICASSO: software for investigating the reliability of ICA estimates by clustering and visualization (Himberg et al., 2004) as implemented in the GIFT software by running ICA several times and selecting the most representative run (Du et al., 2014). ICASSO was used to evaluate component stability. The Infomax ICA algorithm was run for 100 times and clustered together within the ICASSO framework (Himberg et al., 2004). The run with the closest ICs to the centrotypes of stable clusters (ICASSO cluster quality index >0.8) (Iraji et al., 2020) was selected as the best run and used for future analysis (Ma et al., 2011).

We utilized ICA with different model orders (25, 50, 75, and 100) to identify ICNs at multiple spatial scales. ICNs were identified from each model order, and included components with peak activations in gray matter and low-frequency time courses (Calhoun et al., 2009). The subject-specific ICNs time courses were calculated using the spatial multiple regression technique (Calhoun et al., 2004). Before calculating static FNC, an additional cleaning procedure was performed on the time courses of ICNs to reduce the effect of the remaining noise and improve the detection of FNC patterns. The procedures were as follows: (1) remove the linear, quadratic, and cubic trends; (2) regress out the motion realignment variables; (3) replace outliers with the best estimate using a third-order spline fit; and (4) bandpass filter using a fifth-order Butterworth filter with frequency cutoff of 0.01–0.15 Hz (Allen et al., 2011).

FNC analysis at multiple spatial scales

FNC was computed between each pair of ICNs by calculating the Pearson correlation coefficient between ICNs (Allen et al., 2011; Calhoun et al., 2003; Jafri et al., 2008). We calculated a 2D symmetric ICN × ICN FNC matrix for each subject and aggregated the FNC matrix from all subjects into an augmented 2D matrix. We then calculated the mean FNC matrix of all subjects for further analysis.

To investigate the group differences between HC and SZ in FNC, we performed a generalized linear model (GLM). We fit the GLM model with age, gender, data acquisition site, and meanFD as covariates. The meanFD was added to the GLM to account for any residual motion effect that was not removed in the preprocessing steps. The group differences between HC and SZ in FNC were evaluated by the t-value and p-value of statistical comparisons of the GLM.

SVM-based classification

The SVM (Verma and Salour, 1992) is so far the most popular classification method due to its favorable characteristics of high accuracy, ability to deal with high-dimensional data, and versatility in modeling diverse sources of data. The SVM has been widely applied in numerous neuroimaging classification studies and has achieved remarkable results due to its excellent generalization performance. Our motivation of using SVM over other approaches was due to its sensitivity, resilience to overfitting, ability to extract and interpret features, and superior performance in fMRI data classification (De Martino et al., 2008; Ecker et al., 2010; Liu et al., 2013; Pereira et al., 2009; Saha et al., 2021; Vergun et al., 2013; Wang et al., 2019). To investigate the group differences, we applied a binary SVM nonlinear Gaussian RBF kernel classifier (Hsuet al., 2003), due to the fact that the Gaussian RBF kernel generally outperforms other kernels in brain MRI imaging classification (Kumari, 2013; Madheswaran and Anto Sahaya Dhas, 2015). The classification model was trained and crossvalidated using the dataset of 477 HCs and 350 SZs. The figure (Fig. 1) presents the pipeline for training and testing the SVM model.

FIG. 1.

Pipeline of the classification model. The whole dataset was split into a training set (80%) and a testing set (20%). Feature selection was performed on the training set (50% were selected randomly every time). To select the most predictive features, we repeated the feature selection process for ten rounds and retained those features with a high average weight (top 70%) among all the rounds. The final SVM model was built based on the selected FNC features. We ran the modeling process for a total number of 100 iterations to obtain a stable SVM model. FNC, functional network connectivity; SVM, support vector machine. Color images are available online.

Feature selection

Each FNC pair of the functional connectivity among ICNs was considered as the input feature for classification, and the category of group HC or SZ was considered as the response vector. Given that some of those FNC features might be noninformative or redundant for classification, we performed feature selection using Relief (Verma and Salour, 1992) to improve classification performance and to speed up computation. Relief ranks predictors with k nearest neighbors (we set k to 10 for calculation simplicity). The function returns the indices of the most important predictors and the weights of the predictors. Feature selection was carried out before classifier training through the recursive feature elimination step. For each round of feature selection, 50% of the training data (as shown in Fig. 1) was selected. We repeated it for 10 rounds and retained those features with a high average weight (top 70%) among all the rounds. We thus narrowed the set of features to a subset of the original feature set, which eliminated the noninformative and redundant FNC features. The SVM model was built based on the selected FNC feature set.

Recursive validation

To obtain a stable performance of the SVM model, we applied recursive validation. For each iteration, we randomly split the whole dataset into 80% of the training set and 20% of the testing set. The test set was held out for final evaluation. We ran the modeling process for a total number of 100 iterations and evaluated the SVM model based on average specificity, sensitivity, and F1 score (the harmonic mean of the precision and recall) across all iterations.

Results

ICA analysis

Spatial maps of selected ICs at different model orders are shown in the figure (Fig. 2a–d). A total of 127 ICNs were determined between all model orders (15, 28, 36, and 48 from 25, 50, 75, and 100, respectively). ICNs were grouped into functional domains, including cerebellum, cognitive control, default mode, somatomotor, subcortical, temporal, and visual (Fig. 3).

FIG. 2.

ICN maps selected from model order of 25 (a), 50 (b), 75 (c), and 100 (d). A total number of 127 ICNs were determined between all model orders (15, 28, 36, and 48 from 25, 50, 75, and 100, respectively). All the ICNs were identified from each model order and included components with peak activations in gray matter as well as low-frequency time courses. ICN, intrinsic connectivity network. Color images are available online.

FIG. 3.

Identified ICNs across all model orders 25, 50, 75, and 100. ICNs were divided into groups (functional domains) based on their anatomical and functional properties, and include CR, CC, DM, SM, SB, TP, and VS. Each functional domain is displayed at the three most informative slices. CC, cognitive control; CR, cerebellum; DM, default mode; SB, subcortical; SM, somatomotor; TP, temporal; VS, visual. Color images are available online.

The results show that at lower model orders (25–50), signal sources tend to merge into singular components, which then split into several subcomponents at higher model orders. These findings are consistent with previous and most recent findings (Abou-Elseoud et al., 2010; Rachakonda et al., 2017). Figures (Fig. 2c, d) show that higher model order ICA (75–100) tends to parcellate the brain into focal functionally homogeneous distinct regions, which are consistent with previous studies (Calhoun and Adali, 2012; Calhoun and de Lacy, 2017; Iraji et al., 2009). The ICN templates (as shown in Figs. 2 and 3) are openly shared at http://trendscenter.org/data.

FNC analysis

The figures (Fig. 4 a, b) present the mean FNC (z-fisher score) between 127 ICNs between all model orders. The figure (Fig. 4a) was sorted by ICNs, and then model orders within each domain (called an FNC “block plot”), and the figure (Fig. 4b) was sorted by model order and then the domain (called an FNC “finger plot”), aiming to visualize patterns existing at and between different model orders. In the block plot, we observed that the cerebellum domain was highly correlated with itself, anticorrelated across all but subcortical domains. The somatomotor, subcortical, visual (and to a lesser degree temporal) domains were more homogeneous than the default mode, cerebellum, and cognitive control domains both within and between model orders. We also observed strong anticorrelations between default mode versus somatomotor and temporal between model orders 25 and 100 in the finger plot, which were not clearly seen in the other model orders.

FIG. 4.

(a) (Left). Average FNC plots. We calculated the mean FNC (z-fisher score) based on the aggregated FNC matrix of all individuals. (a) (Left). Block plot of mean FNC matrix between model orders. The ICNs in this FNC matrix were sorted by domains first, and within each domain, ICNs were sorted by model orders (from 25 to 100). The dotted lines in each domain divide different model orders (b) (Right). Finger plot of mean FNC matrix between model orders. ICNs within each model order were sorted in the order of CR, CC, DM, SM, SB, TP, and VS. As we can see in the plot, DM from model order 25 shows strong anticorrelations with SM and TP of model order 100; similar anticorrelations were also observed between model orders 75 and 100, which did not show up in the other sections of the matrix. Color images are available online.

We also evaluated the group differences between SZ and controls (Fig. 5a, b). Strong increases were observed between the ICNs of the subcortical domain and the ICNs of the temporal, visual, and somatomotor domains in the SZ group. A similar increased FNC in SZ was observed between the ICNs of the cerebellum domain and the ICNs of somatomotor, temporal, and visual domains. The most dominant decreases in FNC in the SZ group were found mainly within the ICNs of the somatomotor domain, and between the ICNs of the somatomotor domain and the ICNs of the temporal and visual domains. Furthermore, slight decreases in SZ were seen between the cerebellum and subcortical compared with the HC group.

FIG. 5.

(a) (Left). Block plot of GLM contrast of difference mean FNC matrix (HC–SZ). To test the group difference between HC and SZ, we fit a GLM with age, gender, data acquisition site, and meanFD as covariates. This figure shows the difference between HC and SZ in FNC. Blue areas indicate increased FNC in SZ compared with HC, and red areas indicate decreased FNC in SZ compared with HC. The statistical results of p-values of the GLM were corrected for multiple comparisons using a 5% FDR (b) (Right). Finger plot of GLM contrast of difference mean FNC matrix (HC–SZ). (a, b) Show the intensity [−sign(T) × log10(FDR)] matrixes of both block and finger plots for the mean FNC of the GLM model, where T is the t-statistic values of the GLM. The lower triangles show covariate pairs that were significantly different (FDR <0.05). FD, framewise displacement; FDR, false discovery rate; GLM, generalized linear model; HC, healthy control. Color images are available online.

SVM model

Figure 6a and b represents the block plot and the finger plot of average weights of FNC features of 100 iterations. This gives us a better view of the importance of each FNC feature in detecting the group differences. As shown in the figures (Fig. 6a, b), FNC features in the cerebellum, somatomotor, subcortical, temporal, and visual domains consistently contribute to the classification results for both within and between different model orders. These findings suggest highly distinct functional roles for these ICNs in differentiating between HCs and SZs. It is noticeable that the FNC within the somatomotor and visual has contributed most to the classification. The predictive abilities of the FNC features within those two domains were evenly distributed at all model orders, as they are shown in the diagonal regions within somatomotor and visual (in Fig. 6a). Generally, strong predictive abilities of FNC features were observed in the FNC pairs between the subcortical domain and the domains of the cerebellum, the somatomotor, and the temporal across all model orders. Similarly, FNC features between the somatomotor and the domains in the cerebellum and temporal were generally predictive.

FIG. 6.

(a) (Left). SVM block plot of FNC feature average weights of 100 iterations. The hotspots indicate FNC features with strong weights in predicting group differences of HC and SZ. Statistically, the relevance level of a relevant feature is expected to be >0 and that of an irrelevant one is expected to be 0 (or negative). The cerebellum, somatomotor, subcortical, and temporal domains almost always contribute to the classification. Generally, strong predictive abilities of FNC features in the somatomotor and visual domains were seen at all model orders. Visual was predictive only within the region between somatomotor and visual, but only seen at high model order (100) (b) (Right). SVM finger plot of FNC feature average weights of 100 iterations. Color images are available online.

To better understand the between-model order differences, we computed a two-sample t-test (Cressie and Whitford, 1986), to identify significant differences within and between model orders. For each iteration, we calculated an average feature weight for each model order; we then aggregated all the average feature weights across 100 iterations. The two-sample t-tests were performed between all pairs of model orders using their aggregated feature weights. It (Fig. 7) shows the statistical comparison (intensity) plot between model orders. A negative intensity value in the lower triangle of the figure indicates that the average feature weights of model orders in the row were smaller than the ones in the column, and a positive intensity value indicates that the model orders' average feature weights in the row were larger than the ones in the column. As shown in the figure (Fig. 7), features in the between model orders of 25–50 were consistently more predictive than features in other within and between model orders. Between-model order features were consistently more predictive than higher individual model orders (75–100), but less predictive than lower model orders (25–50). And we also see that features in the model order of 50 are more predictive than most other model orders, except the between model orders 25–50. Features in the between model orders 75–100 are consistently less predictive than most within and between model orders, except model order of 100.

FIG. 7.

The intensity [−sign(T) × log₁₀(FDR)] of average feature weights between model orders. We computed a two-sample t-test between the averaged feature weights of every two model orders, to identify significant differences between them. The negative intensity value in the lower triangle of the figure indicates the average feature weights of model orders in the row (x) were smaller than the ones in the column (y), and the positive intensity value in the lower triangle indicates the model orders' average feature weights in the row were larger than the ones in the column. The upper triangle shows significant differences (p < 0.05, FDR corrected) between each pair of model orders, whose direction is opposite to the lower triangle. Color images are available online.

To highlight the variations between domains of model orders, we also calculated the average feature weights (Fig. 8a) and maximum feature weights (Fig. 8b) in the domain levels. Within each domain, feature weights were aggregated by model orders from 25 to 100. It shows that higher average feature weights were mainly seen in the somatomotor, subcortical, temporal, and visual domains. Strong predictive abilities of FNC features within the somatomotor and visual (highlighted in the figures) were seen at all model orders.

FIG. 8.

(a) (Left). Average feature weight between domains. Each domain contains four model orders, from 25, 50, 75 to 100, and all the features within each model order were averaged across 100 iterations (b) (Right). Maximum feature weight within each domain. Each domain column contains four averaged features of different model orders, from 25, 50, 75 to 100. The maximum feature weight within each domain was selected based on the averaged feature weight across 100 iterations. It shows that higher average feature weights are mainly seen in somatomotor, subcortical, temporal, and visual. Color images are available online.

We then performed the two-sample t-test between domains of different model orders (Cressie and Whitford, 1986), to identify significant differences between domains. We calculated the intensity values between domains of different model orders (as shown in Fig. 9). For each ICN in the figure, we averaged all the feature weights within each model order of 25, 50, 75, and 100. It indicates that features in somatomotor were always predictive than features in other domains at all model orders, such as cerebellum, cognitive control, and default mode (blue regions in the somatomotor row of the lower triangle), and subcortical, temporal, and visual (hot regions in the somatomotor column of the lower triangle). Besides, features in visual were generally predictive than other domains at all model orders, except somatomotor. Features of lower model orders (25–50) in subcortical were generally less predictive than features in other domains.

FIG. 9.

The intensity [−sign(T) × log₁₀(FDR)] of average feature weights between domains. The negative intensity values in the lower triangle indicate smaller average feature weights of the ICN in the row (x) compared with the ones in the column (y), and the positive intensity values in the lower triangle indicate larger average feature weights of ICN in the row compared with the column. The upper triangle shows the significant differences (p < 0.05, FDR corrected) between each pair of domains, whose direction is opposite to the lower triangle. Color images are available online.

We also evaluated the performance of the SVM model. As shown in Table 1, the average accuracy of the SVM model using 100 iterations was 77.4% (outperformed the null accuracy 57.6%), with a specificity of 85.2%, 67.2% sensitivity, and 71.5% F1 score.

Table 1.

The Average Performance of Support Vector Machine Model for 100 Iterations

	Accuracy	Specificity	Sensitivity	F1
Support vector machine	0.7743	0.8522	0.6717	0.7153

Finally, we compared the predictive accuracies of SVM models built by different model orders. We individually trained the SVM model using each number of model order, and then we computed a two-sample t-test between the accuracies of every two SVM models of different model orders. As shown in the figure (Fig. 10), higher individual model orders always have higher accuracies compared with lower individual model orders. This is due to the fact that including more model orders always resulted in the increase of predictive features. It is also observed that the between model orders of 50–75, 50–100, and 75–100 outperformed individual model orders of 25, 50, 75, and 100, while individual model orders (50–100) outperformed lower between model orders of 25–50 and 25–75.

FIG. 10.

The intensity [−sign(T) × log₁₀(FDR)] of SVM accuracies between model orders. The negative intensity value in the lower triangle of the figure indicates the average feature weights of model orders in the row (x) were smaller than the ones in the column (y), and the positive intensity value in the lower triangle indicates the model orders' average feature weights in the row were larger than the ones in the column. The upper triangle shows the significant differences (p < 0.05, FDR corrected) between each pair of model orders, whose direction is opposite to the lower triangle. The top circle indicates that higher individual model orders always have higher accuracies compared with lower individual model orders. The middle circle shows that model orders of 50–100 outperform lower between model orders. The bottom circle indicates that higher between model orders outperform individual model orders. Color images are available online.

Discussion

In this article, we introduce a multimodel order ICA approach to estimate FNCs at multispatial scales. We investigated whether multispatial-scale FNC features were able to discriminate between SZ individuals and HC groups, and found both consistency and uniqueness of FNC patterns within and between model orders. Importantly, we revealed that additional information can be preserved in the between-order FNC that might be ignored in single-model order analysis. In addition, some interesting findings were only observed at lower model orders or higher model orders. Specifically, the HC–SZ differences were observed between subcortical versus temporal and subcortical versus cerebellum mostly at lower model order FNC (25–75), but between somatomotor versus visual networks were mainly at higher model order (100). Results highlighted findings we learned from multiple-model orders that would have been missed otherwise.

In addition, we see consistency between the identified significant group differences in FNC and the SVM feature weights (as shown in Figs. 5 and 6). We observed that HC generally showed high FNC in somatomotor, visual at all model orders compared with SZ (Fig. 5a), and these cells were also found to be the most predictive features in the SVM model (Fig. 6a). For between model orders, we observed that although somatomotor versus subcortical and visual versus subcortical were relatively high in SZ compared with HC, they were less predictive than somatomotor and visual. It is also interesting to note that some noticeable differences between HC and SZ in the difference mean FNC matrix (Fig. 5a) do not necessarily lead to their predictive abilities (the same areas in Fig. 6a are not highlighted). For example, SZs were generally higher than HCs in cerebellum versus visual and cerebellum versus temporal in the difference mean FNC matrix (Fig. 5a), but they were not noticeably predictive in the SVM model (Fig. 6a).

Some other useful findings were also indicated by jointly studying the functional brain imaging data within and between the different number of model orders. By comparing the FNC feature predictive ability between different model orders, we surprisingly found that the between model orders always maintained useful information in detecting group differences, particularly the features between model orders of 25–50 were the most SZ discriminative ones compared with features in other within and between model orders. And in many cases, features in the between model orders were more predictive than individual model order according to our observations, suggesting that investigation of between-model order FNC may provide informative details that could be missed otherwise. Results also highlighted the differentiating power of somatomotor and visual domains on classifying HC versus SZ both within and between model orders, compared with other domains.

There are several limitations to this study. First, to obtain a stable feature set, we performed 10 rounds of feature selection within each iteration when turning the SVM model, and within each round of feature selection, we implemented a strategy similar to bootstrap sampling (Breiman, 1996), by selecting a subset (50%) of the training data for feature selection. In future work, the performance of the classification model can be improved by increasing the rounds of feature selection and by exploring a range of selected subsets for feature selection, to get a more stable feature set. Second, we selected the top 70% ranked FNC features to build the SVM model, based on our previous experience, that in general, selecting top 70% ranked features maintains the maximum predictive information of the dataset by eliminating redundant features. Going forward, we can improve the performance of the classification model by exploring a range of selected features, to get an optimal predictive model. We tried both using feature selection and not using it to build the SVM model; the results showed that the performance of using feature selection was slightly better than not using it. In addition to the SVM model, we tried another classification model random forest (Ho, 1995); the performance of SVM was slightly better compared with the random forest for our case, so we keep SVM results. For future work, more efforts will be put in to explore other machine learning techniques, such as boosting (Schapire, 2003), to improve the accuracy of the learning classifier.

Conclusions

In summary, based on fMRI data from 477 healthy participants and 350 SZ patients, we compared the group differences obtained at different model orders by evaluating the multiorder FNCs, as well as testing the classification power of FNC at different model orders. A comprehensive visualization of the relationships through FNC block plots and FNC finger plots was introduced. Results suggest that additional information was observed both at different model orders and in between-model order relationships, which goes beyond the known general modular relationships in ICNs. We are also releasing four model order templates to the public for future use. Our study highlights the benefit of studying FNC within and between multiple spatial scales. This work expands upon and adds to the relatively new literature on resting fMRI-based classification and prediction.

Footnotes

Acknowledgments

We thank the authors of the studies for collection of the original data.

Authors' Contributions

X.M. contributed to conceptualization; formal analysis; investigation; methodology; visualization; writing—original draft; writing—review and editing. A.I. performed conceptualization; investigation; methodology; visualization; writing—review and editing. Z.F. performed data collection and preprocessing. P.K., A.B., J.F., S.M., D.H.M., B.A.M., G.P., S.G.P., A.P., J.T., and T.v.E. contributed to data collection and resources. J.S. performed writing—review and editing. V.D.C. contributed to conceptualization; funding acquisition; investigation; methodology; project administration; resources; writing—review and editing.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

National Institutes of Health Awards R01MH118695, R01MH117107, and RF1AG063153 were honored to V.D.C.

References

Abou-Elseoud

, Starck

, Remes

, et al. 2010. The effect of model order selection in group PICA. Hum Brain Mapp, 31:1207–1216.

Allen

, Erhardt

, Damaraju

, et al. 2011. A baseline for the multivariate comparison of resting-state networks. Front Syst Neurosci, 5:1–23.

Anderson

, Cohen

. 2013. Decreased small-world functional network connectivity and clustering across resting state networks in Schizophrenia: an fMRI classification tutorial. Front Hum Neurosci, 7:1–18.

Beckmann

, DeLuca

, Devlin

, et al. 2005. Investigations into resting-state connectivity using independent component analysis. Phil Trans R Soc B Biol Sci, 360:1001–1013.

Breiman

1996. Bagging predictors. Mach Learn, 24:123–140.

Calhoun

, Adali

. 2012. Multisubject independent component analysis of fMRI: a decade of intrinsic networks, default mode, and neurodiagnostic discovery. IEEE Rev Biomed Eng, 5:60–73.

Calhoun

, Adali

, Hansen

, et al. 2003. ICA of functional MRI data: an overview, no. April:. 281–288. [Epub ahead of print]; DOI: 10.1.1.3.7473.

Calhoun

, Adali

, Pearlson

, et al. 2001. A method for making group inferences from functional MRI data using independent component analysis. Hum Brain Mapp, 14:140–151.

Calhoun

, de Lacy

. 2017. Ten key observations on the analysis of resting-state functional MR imaging data using independent component analysis. Neuroimaging Clin N Am, 27:561–579.

10.

Calhoun

, Eichele

, Pearlson

. 2009. Functional brain networks in Schizophrenia: a review. Front Hum Neurosci, 3:1–12.

11.

Calhoun

, Kiehl

, Pearlson

. 2008. Modulation of temporally coherent brain networks estimated using ICA at rest and during cognitive tasks. Hum Brain Mapp, 29:828–838.

12.

Calhoun

, Pekar

, Pearlson

. 2004. Alcohol intoxication effects on simulated driving: exploring alcohol-dose effects on brain activation using functional MRI. Neuropsychopharmacology, 29:2097–2107.

13.

Cressie

NAC

, Whitford

. 1986. How to use the two sample T-test. Biometrical J, 28:131–148.

14.

Damoiseaux

, Beckmann

, Sanz Arigita

, et al. 2008. Reduced resting-state brain activity in the ‘Default Network’ in normal aging. Cereb Cortex, 18:1856–1864.

15.

De Martino

, Valente

, Staeren

, et al. 2008. Combining multivariate voxel selection and support vector machines for mapping and classification of fMRI spatial patterns. Neuroimage, 43:44–58.

16.

, Ma

, Fu

, et al. A novel approach for assessing reliability of ICA for FMRI analysis. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Florence, Italy, 2014, pp. 2084–2088.

17.

, Fu

, Sui

, et al. 2020. NeuroMark: an automated and adaptive ICA based pipeline to identify reproducible fMRI markers of brain disorders. Neuroimage, 28:102375.

18.

Ecker

, Rocha-Rego

, Johnston

, et al. 2010. Investigating the predictive value of whole-brain structural MR scans in autism: a pattern classification approach. Neuroimage, 49:44–56.

19.

Himberg

, Hyvärinen

, Esposito

. 2004. Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage, 22:1214–1222.

20.

TK.

1995. Random decision forests. In Proceedings of the International Conference on Document Analysis and Recognition, ICDAR, vol. 1. Montreal, QC: IEEE, pp. 278–282.

21.

Hsu

C-W

, Chang

C-C

, Lin

C-J

. 2003. A Practical Guide to Support Vector Classification. 1396–1400. [Epub ahead of print]; DOI: 10.1177/02632760022050997.

22.

Iraji

, Faghiri

, Fu

, et al. 2021. Multi-spatial scale dynamic interactions between functional sources revel sex-specific changes in Schizophrenia. Netw Neurosci. 1–48. [Epub ahead of print]; DOI: 10.1162/netn_a_00196.

23.

Iraji

, Faghiri

, Lewis

, et al. 2020. Tools of the trade: estimating time-varying connectivity patterns from fMRI data. Social Cogn Affect Neurosci. 1–26. [Epub ahead of print]; DOI: 10.1093/scan/nsaa114.

24.

Iraji

, Faghiri

, Lewis

, et al. 2009. Ultra-High-Order ICA: An Exploration of Highly Resolved Data-Driven Representation of Intrinsic Connectivity Networks (Sparse ICNs), vol. 11138. International Society for Optics and Photonics. [Epub ahead of print]; DOI: 10.1117/12.2530106.

25.

Iraji

, Fu

, Damaraju

, et al. 2019. Spatial dynamics within and between brain functional domains: a hierarchical approach to study time-varying brain function. Hum Brain Mapp, 40:1969–1986.

26.

Jafri

, Pearlson

, Stevens

, et al. 2008. A method for functional network connectivity among spatially independent resting-state components in Schizophrenia. Neuroimage, 39:1666–1681.

27.

Kumari

2013. SVM classification an approach on detecting abnormality in brain MRI images. Int J Eng Res Appl (IJERA), 3:1686–1690.

28.

Lewis

, Bosman

, Fries

. 2015. Recording of brain activity across spatial scales. Curr Opin Neurobiol, 32:68–77.

29.

Liu

, Tosun

, Weiner

, et al. 2013. Locally linear embedding (LLE) for MRI based Alzheimer's disease classification. Neuroimage, 83:148–157.

30.

, Correa

, Li

, et al. 2011. Automatic identification of functional clusters in fMRI data using spatial dependence. IEEE Trans Biomed Eng 58 (12 PART 1):3406–3417.

31.

Madheswaran

, Anto Sahaya Dhas

. 2015. Classification of brain MRI images using support vector machine with various kernels. Biomed Res (India), 26:505–513.

32.

Pariyadath

, Stein

, Ross

. 2014. Machine learning classification of resting state functional connectivity predicts smoking status. Front Hum Neurosci 8. [Epub ahead of print]; DOI: 10.3389/fnhum.2014.00425.

33.

Perdue

, Diamond

. 2013. Effects of spatial pattern scale of brain activity on the sensitivity of DOT, fMRI, EEG and MEG. PLoS One 8. [Epub ahead of print]; DOI: 10.1371/journal.pone.0083299.

34.

Pereira

, Mitchell T Botvinick

. 2009. Machine learning classifiers and fMRI: a tutorial overview. Neuroimage 45 (1 Suppl):S199–S209.

35.

Rachakonda

, Du

, Calhoun

. 2017. Model order prediction in ICA. In Poster in Organization for Human Brain Mapping (OHBM). Vancouver, Canada.

36.

Saha

, Damaraju

, Rashid

, et al. 2021. A classification-based approach to estimate the number of resting functional magnetic resonance imaging dynamic functional connectivity states. Brain Connect, 11:132–145.

37.

Salvador

, Sarró

, Gomar

, et al. 2010. Overall brain connectivity maps show cortico-subcortical abnormalities in Schizophrenia. Hum Brain Mapp, 31:2003–2014.

38.

Schapire

RE.

2003. The boosting approach to machine learning: an overview BT—nonlinear estimation and classification. In: Denison DD, Hansen MH, Holmes CC, Mallick B, Yu B. New York, NY: Springer New York, pp. 149–171. [Epub ahead of print]; DOI: 10.1007/978-0-387-21579-2_9.

39.

Vergun

, Deshpande

, Meier

, et al. 2013. Characterizing functional connectivity differences in aging adults using machine learning on resting state fMRI data. Front Comput Neurosci, 7:1–20.

40.

Verma

, Salour

. 1992. The feature selection problem: traditional methods and a new algorithm. AAAI-92 Proc, 256:175–200.

41.

Wang

, Li

, Zhang

, et al. 2019. Support vector machine for analyzing contributions of brain regions during task-state fMRI. Front Neuroinform, 13:1–12.

42.

Whitfield-Gabrieli

, Thermenos

, Milanovic

, et al. 2009. Hyperactivity and hyperconnectivity of the default network in Schizophrenia and in first-degree relatives of persons with Schizophrenia. Proc Natl Acad Sci U S A, 106:1279–1284.

43.

, Du

, Chen

, et al. Information Technology, New Haven, and New Haven. 2018. Comparing brain graphs in which nodes are regions of interest or independent components: a simulation study. J Neurosci Methods. 61–68. [Epub ahead of print]; DOI: 10.1016/j.jneumeth.2017.08.007.Comparing.