A telescopic independent component analysis on functional magnetic resonance imaging dataset

A macroscopic brain’s function can be conceptualized as coordination and interaction among functional sources, a model incorporating the fundamental principles of segregation and interaction (Genon, Reid, Langner, Amunts, & Eickhoff, 2018; Iraji et al., 2022). From this perspective, the brain is composed of distinct functional sources that dynamically interact with each other. These functional sources can be characterized by temporally synchronized patterns, which can be assessed at the macroscale using resting-state functional magnetic resonance imaging (rs-fMRI) and techniques such as independent component analysis (ICA; Calhoun & Adali, 2012; van den Heuvel & Hulshoff Pol, 2010). ICA is a robust blind source separation technique that separates complex rs-fMRI signals into spatial patterns that are temporally coherent and maximally independent of each other (Calhoun & Adali, 2012; Calhoun, Adali, Pearlson, & Pekar, 2001; Calhoun & de Lacy, 2017). These spatial patterns and associated time courses formed the estimation of functional sources and are commonly referred to as intrinsic connectivity networks (ICNs; Allen et al., 2014). Importantly, these ICNs can be captured at different spatial scales (Abou-Elseoud et al., 2010; Iraji et al., 2022, 2023). Previous research has underscored the significance of the analysis of ICNs at multiple spatial scales (Iraji et al., 2022; Li, Zhu, & Fan, 2018; Meng et al., 2022, 2023). Previous studies that explore multiple spatial scales can be divided into two major categories. The category relies on a singular set of nodes, such as predefined regions or a single-model-order ICA, to reconstruct different scales of functional hierarchy using various modularity or clustering techniques (Doucet et al., 2011); however, these approaches have no information that can be explained by looking at the other spatial scales. The other category of approaches appreciates that ICNs at different spatial scales can contain distinctive functional sources with their unique information. For instance, previous studies utilized ICA with different model orders to capture ICNs at different spatial scales to study whole-brain functional connectivity (Iraji et al., 2022; Meng et al., 2022). Lower model orders yield large-scale, spatially distributed ICAs (Beckmann, DeLuca, Devlin, & Smith, 2005; Calhoun, Kiehl, & Pearlson, 2008; Iraji et al., 2016), whereas higher model orders result in more finely grained, spatially specific ICNs (Abou-Elseoud et al., 2010; Allen et al., 2011). Importantly, the interaction between large-grained scale and small-grained scale can convey crucial insights into the brain’s functioning. For instance, multiscale brain analysis better captures sex-specific schizophrenia changes (Iraji et al., 2022). However, these techniques require post hoc analysis to construct the spatial hierarchy and assign ICNs to different spatial scales, as the method does not impose direct constraints on the size of ICNs.

Here, we present a new ICA strategy, called telescopic ICA (TICA), to construct spatial functional hierarchy and estimate ICNs across multiple spatial scales. The approach is a recursive technique that leverages ICNs from the previous scale to guide the ICA to obtain ICNs for another scale. Notably, our approach is designed to show how a large-scale network decomposes through a small-scale network using zooming in based on the prior level of the functional source.

Schizophrenia represents a psychotic disorder characterized by different cognitive impairments and a decline in both social and occupational functioning (American Psychiatric Association, 2013; van den Heuvel & Fornito, 2014). This condition is multifaceted, is diagnosed through a syndromic approach, and involves excluding other potential diagnoses. Schizophrenia lacks distinctive symptoms and is primarily diagnosed based on clinical observations of positive symptoms (such as delusions, hallucinations, disorganized speech, and disorganized or catatonic behavior) as well as negative symptoms (including apathy, blunted affect, and anhedonia), along with a noticeable deterioration in social functioning (American Psychiatric Association, 2013). Schizophrenia exhibits a substantial overlap with both schizoaffective disorder and psychotic bipolar disorder, in both symptoms, genetic factors, and other biomarkers (Clementz et al., 2016). Schizophrenia is theorized to be a developmental disorder marked by disrupted brain function, characterized by either functional dysconnectivity or alterations in functional integration (Friston & Frith, 1995; Kahn et al., 2015). Research also has elucidated sex differences in the incidence and clinical manifestation of mental disorders (Dalsgaard et al., 2020; Paus, Keshavan, & Giedd, 2008). Females with schizophrenia tend to show more depressive symptoms, whereas males often experience more negative symptoms (Ochoa, Usall, Cobo, Labad, & Kulkarni, 2012). Also, recent studies highlight the importance of analyzing the brain in multiple scales to capture sex-specific schizophrenia changes (Iraji et al., 2022), and previous studies have shown high-order functional connectivity applied to the diagnosis of psychotic disorders (Herzog et al., 2022; Li et al., 2023). Therefore, studying brain functional sources at different scales can provide crucial information about brain functional integration and its schizophrenia changes, potentially improving our understanding of the actual brain pathology underlying different schizophrenia subcategories. Here, we used the TICA approach to identify the spatial map of the default mode network (DMN) in two scales and provide important information about brain functional sources and its schizophrenia changes, potentially improving our understanding of the actual brain pathology and, thus, eventually improving treatments and care for individuals with schizophrenia. Our results show that the TICA approach can capture significant group differences in the DMN area between healthy people (HC) and schizophrenia patients (SZ) that are missed in single-scale ICA.

We utilized a dataset from the Functional Imaging Biomedical Informatics Research Network (FBIRN) for our investigation (Potkin & Ford, 2009). The FBIRN dataset was collected from seven sites. The same rs-fMRI parameters were used across all sites: a standard gradient echo-planar imaging (EPI) sequence, repetition time (TR)/echo time = 2,000/30 ms, voxel spacing size = 3.4375 × 3.4375 × 4 mm, slice gap = 1 mm, flip angle = 77°, field of view = 220 × 220 mm, and a total of 162 volume. Six of the seven sites used 3-T Siemens Tim Trio scanners, and one used a 3-T General Electric Discovery MR750 scanner. The dataset was chosen using the following quality control: (a) individuals diagnosed with either typical control or schizophrenia conditions, (b) data demonstrating high-quality registration to the EPI template, (c) minimal head motion characterized by less than 3° rotations and 3-mm translations in every direction, and (d) mean framewise displacement of less than 0.2. Consequently, our analysis was performed on rs-fMRI data derived from 310 individuals, comprising 150 HC and 160 SZ, as summarized in Table 1. All participants are also at least 18 years old and gave written informed consent before enrollment.

The rs-fMRI data preprocessing was conducted using the Statistical Parametric Mapping toolbox (SPM12; https://www.fil.ion.ucl.ac.uk/spm/). Preprocessing steps include discarding the initial five volumes, rigid motion correction, and slice-time correction. Subsequently, individual subject data were registered to an MNI EPI template, resampled, and smoothed with a 6-mm Gaussian kernel, and variance normalization was applied on voxel time (Iraji et al., 2022).

We used the Group Independent Component Analysis Toolbox (GIFT) toolbox (https://trendscenter.org/software/gift/; Calhoun & Adali, 2012; Calhoun et al., 2001; Iraji et al., 2021) to apply the TICA. To begin, the 30 PCs that explained the maximum variance of each subject and the 20 group-level PCs that explained the maximum variance of each estimator-specific dataset were used as the input for spatial GICA. We selected a spatial GICA model order of 20 to obtain large-scale ICNs (Iraji et al., 2016, 2019). For the GICA algorithm, we used the infomax algorithm (Bell & Sejnowski, 1995; Lee, Girolami, & Sejnowski, 1999), ran it 100 times, and utilized the ICASSO framework to identify the best estimate, defined as the independent component closest to the stable cluster (Himberg, Hyvärinen, & Esposito, 2004). We employed back-reconstruction using a GICA technique (Calhoun et al., 2001) to compute subject-specific independent components and time courses because it provides more accurate time courses and spatial maps.

Afterward, we employed a selection process for identifying ICNs. We identified ICNs using prior knowledge and the following criteria: (a) dominant low-frequency fluctuations in their time courses, determined through dynamic range assessment and the low-frequency to high-frequency power ratio; (b) peak activation located within gray matter regions; (c) minimal spatial overlap with vascular and ventricular structures; and (d) low spatial similarity with motion and other recognized artifacts (Allen et al., 2011). Next, we picked the DMN, VN, and RFPN and performed the weighting step for all subjects’ BOLD signals with their spatial map of Scale 1. Finally, we applied the GICA to the weighted BOLD signals. Similar to the previous scale, we applied two levels of PCA. Six subject-level and four group-level PCs that explained the maximum variance of the dataset were used to input for the GICA. Previous works (Abou-Elseoud et al., 2010; Allen et al., 2011; Hu et al., 2016; Iraji et al., 2022; Li, Adali, & Calhoun, 2007) highlight the impact of model order on ICA results; we used elbow point criteria to determine the optimal number of PCs. All other parameters in the GICA to extract ICN Scale 2 are the same as those for Scale 1. We also applied GICA Scale 2 with model orders 2, 3, 5, and 6 on the DMN to investigate how model-order choice affects the results of the TICA.

We evaluated the replicability (Adali et al., 2022; National Academies of Sciences, Engineering, and Medicine, 2019) of our approach using the initial random division of the dataset into two distinct subsets, each of which underwent independent processing through our pipeline. This division and processing were repeated 100 times, with each iteration yielding one DMN Scale 1 and its corresponding four zoomed-in DMNs extracted from Scale 2. To quantify the replicability of our approach, we calculated the average of spatial similarity between each pair of ICNs from two independent half splits using Pearson correlation (Figure 2).

We also assessed the similarity between the individual-level DMN and the group-level DMN in Scales 1 and 2 as well as the similarity between the individual-level DMN in Scale 1 and that in Scale 2. To accomplish this, we obtained the DMN for each participant, and we computed the Pearson correlation between the group-level DMN and the corresponding individual-level DMN to get the similarity at the group and individual levels. We also calculated the Pearson correlation between individual subjects and obtained the similarity at the individual level. This comparative analysis was performed for both Scales 1 and 2 of the DMN.

Moreover, we conducted an additional analysis to assess the test-retest reliability of the TICA to specifically measure individual differences. We used the Human Connectome Project (HCP; Van Essen et al., 2013) datasets. The HCP dataset consists of 1,019 subjects and has four complete rs-fMRI scan sessions. Each session has 1,200 time points, a TR of 0.72 s, and an original voxel size of 2.0 mm isotropic. We ran the TICA separately on two different scan sessions and calculated the Pearson correlation similarity between group-level DMNs in Scales 1 and 2. We also calculated the Pearson correlation similarity for each subject between two different scan sessions.

To assess ICN differences between HC and SZ, we first identified the voxels with a strong contribution to a given ICA (Z > 1.96). We next conducted voxelwise two-sample t tests and applied an adaptive thresholding technique based on the Gamma-Gaussian mixture model (Gorgolewski, Storkey, Bastin, & Pernet, 2012) to identify voxels with significant diagnostic effects.

Moreover, to evaluate the consistency of group differences, we calculated the group difference using a random 80% subset of the dataset. We then measured the similarity between the T-map generated from the entire dataset and the T-map from this 80% subset. This process was repeated 50 times.

To investigate the comparison of DMN components derived through the TICA method using a single-scale ICA, we conducted an ICA with a high model order on the FBIRN dataset. The GICA steps are as follows. First, 120 PCs that explained the maximum variance of each subject were retained for further analysis. Next, 80 group-level PCs that explained the maximum variance of each estimator-specific dataset were used as the input for the GICA. We applied model order 80 as GICA high order. The model order choice was motivated by the TICA procedure, where we initially employed the GICA with a model order of 20 to get networks in Scale 1; then, an additional GICA with a model order of 4 was used to get smaller networks in Scale 2. For the high-model-order GICA algorithm, similar to the TICA approach, we used the infomax algorithm, the ICASSO framework, and GICA back-reconstruction techniques to compute subject-specific independent components and time courses. Afterward, we employed a selection process for identifying the DMN like the TICA approach.

To quantify comparing the results between the TICA and the high-order GICA, we calculated two different approaches. First, we evaluated the similarity between spatial maps generated by the high-order ICA and TICA approaches using Pearson correlation. Second, we conducted a paired-samples t test on spatial maps to compare the TICA approach and the high-model-order ICA across DMNs.

We also applied group comparison analysis for the DMN derived from the high-order GICA. Similar to the TICA approach, we identified the group effect in voxels that significantly contribute to the DMN in HC and SZ groups using an adaptive thresholding technique to identify regions with significant diagnostic effects. Moreover, to evaluate the consistency of group differences, we calculated the group difference using a random 80% subset of the dataset. We then measured the similarity between the T-map generated from the entire dataset and the T-map from this 80% subset. This process was repeated 50 times to find similarity between the T-maps for all components.

According to the criteria specified in the Materials and Methods section, we identified the DMN, VN, and RFPN extracted from Scale 1, along with four zoomed-in ICNs derived from Scale 2. The decomposition of the ICNs across these two spatial scales is visually depicted in Figure 3. In DMN, VN, and RFPN Scale 1, the number of voxels that significantly contribute (Z > 1.96) in spatial maps is 2,948; 3,122; and 3,307, respectively. In Scale 2, DMN significant voxels are 1,077; 2,723; 965; and 1,426, respectively. VN significant voxels are 1,544; 1,239; 2,860; and 1,058, respectively. Significant voxels in RFPN spatial maps are 1,389; 2,777; 968; and 1,130, respectively. These show that the TICA approach captures bigger networks in Scale 1 and smaller networks in Scale 2. In addition, the spatial correlation similarities between DMN Scales 1 and 2 are 0.43, 0.61, 0.42, and 0.61, respectively. The spatial correlation similarities between VN Scales 1 and 2 are 0.72, 0.47, 0.49, and 0.25, respectively. The spatial correlation similarities between RFPN Scales 1 and 2 are 0.27, 0.70, 0.54, and 0.26, respectively. These results display that the telescopic approach can capture distinctive networks in each scale.

We also present the similarity between individual-level DMN spatial maps as well as individual-level and group-level DMN spatial maps generated from the TICA approach using the Pearson correlation coefficient. The computed average similarity between the group and individual levels on Scale 1 was 0.54, and those on Scale 2 were 0.60, 0.72, 0.70, and 0.69, respectively. Also, the average Pearson correlation coefficients at the individual level from DMN Scales 1 and 2 are 0.30, 0.37, 0.52, 0.49, and 0.48. SDs of correlation similarity at the individual level are 0.09, 0.17, 0.18, 0.17, and 0.18, respectively. The visual representation in Figure 4B illustrates the boxplot of these similarities, revealing a higher similarity in DMN Scale 2 as compared with DMN Scale 1.

Figure 5 displays the Pearson correlation coefficient for each subject between two scan sessions evaluating the test-retest reliability of the TICA approach. The results show that the correlation coefficients in Scale 2 are higher than those in Scale 1. In addition, the correlation similarity of group-level DMNs between two scan sessions in Scale 1 is 0.98, and those in Scale 2 are 0.99, 0.99, 0.97, and 0.99, respectively.

To study how model-order choice affects the TICA results, we conducted analyses on TICA Scale 2 with model orders 2, 3, 5, and 6. In general, results (Figure 6) are similar and consistent with the previous ICA work. As model order increases, the components branch out into multiple components (Abou-Elseoud et al., 2010); for instance, the first component in order 2 is formed by combining the first and second components from order 4. The finding also highlights that model orders 3–5 showed spatial overlaps of components.

Figure 7 illustrates the T-map of voxels with significant diagnostic effects (p < 0.01) after adoptive thresholding using the Gamma-Gaussian mixture model. Notably, the maps revealed significant group differences in DMN Scale 2, which were not evident in Scale 1. In DMN Scale 1, the areas exhibiting significant differences include the middle cingulate and paracingulate gyri, posterior cingulate gyrus, cuneus, lingual gyrus, middle occipital gyrus, angular gyrus, and precuneus. Conversely, in DMN Scale 2, significant differences besides the area in Scale 1 are manifested in the superior frontal gyrus, medial orbital, calcarine fissure, surrounding cortex, superior occipital gyrus, and inferior parietal gyrus.

Figure 8 shows consistency in group differences using the Pearson correlation coefficient between T-maps generated from the group difference across the entire dataset and 80% of the dataset over 50 iterations. The results revealed that the average similarity exceeded 0.90 for each component, illustrating the consistency in group differences.

We identified four distinct DMNs extracted through ICA at order 80. Figure 9A displays the spatial maps of these DMNs, denoted as DMN1, DMN2, DMN3, and DMN4. Additionally, Figure 9B presents the outcomes of significant group differences between HC and SZ, derived from all DMNs at a high model order. The identified areas exhibiting significant group differences encompass the inferior parietal gyrus, angular gyrus, precuneus, and cuneus regions. The number of significant voxels in the DMN from ICA order 80 is revealed as 49, 46, 32, and 14, respectively.

We conducted a paired-samples t test to compare TICA and high-model-order approaches across DMNs. Figure 10 displays voxels that show significant differences (p < 0.05) after false discovery rate correction between the two approaches. In addition, we assessed the correlation similarity between the DMN1 high-model-order ICA and the DMN generated from a telescopic approach, revealing values of 0.30, 0.13, 0.35, 0.10, and 0.18. Similar correlation assessments were conducted for the DMN2 high-model-order ICA (0.71, 0.57, 0.15, 0.19, 0.73), DMN3 high-model-order ICA (0.40, 0.01, 0.007, 0.50, 0.44), and DMN4 high-model-order ICA (0.25, 0.10, 0.19, 0.04, 0.08).

Studying brain functional connectivity has improved our understanding of brain functions and the impact of brain disorders. However, current studies in functional connectivity disregard functional connectivity across multiple spatial scales. At best, apply a multimodal ICA (Iraji et al., 2022) to study functional interactions. Still, they do not construct the spatial hierarchy or impose direct constraints on the size of ICNs. In this work, we introduced the TICA as an effective and adaptive tool to identify ICNs, construct spatial functional hierarchies, and estimate functional sources across multiple spatial scales. We also leveraged the approach to study group comparison between HC and SZ, which has been understudied as most research focuses on only a single spatial scale.

Using TICA, we identified five distinct DMNs, including the original network extracted from DMN Scale 1 and its decomposition in four small networks from DMN Scale 2. The ICN Hierarchy section supports the finding that the TICA approach captures bigger networks in the first scale and uses its information to extract smaller networks in the next scale.

We observed that the DMN in the first scale demonstrates a commendable mean reliability of 0.97, highlighting its robustness in effectively detecting the overall DMN. Additionally, the networks in the second scale exhibit reasonably high mean reliabilities, ranging from 0.80 to 0.91, indicating some variability in these more detailed functional sources. Our finding aligns with the observed pattern in a high-model-order ICA approach (Abou-Elseoud et al., 2010).

We observed a more pronounced group effect in Scale 2 in comparison with Scale 1 (see the Clinical Applications section). Additionally, the results indicate that subject-level estimates exhibit greater similarity to both the reference and each other in Scale 2 as opposed to Scale 1 (Figure 4), implying a more accurate estimation of networks. These collective findings suggest that the second scale may be more sensitive in capturing the group effect associated with SZ. Therefore, the TICA approach holds promise in identifying biomarkers for schizophrenia. We observed that DMN3 Scale 2 resembles the spatial maps of the parietal memory network, a network with similar spatial patterns to the DMN but with distinct functional properties (Yang et al., 2014). This finding highlights that the TICA can capture this distinct network and suggests that it might be a network in the DMN hierarchy. Moreover, our analysis identifies significant group differences in the posterior cortex and precuneus areas, echoing similar observations in previous studies (Gilmore, Nelson, & McDermott, 2015; Guo et al., 2020) that contribute to auditory hallucinations in schizophrenia. Additionally, while previous work estimated this network in the first layer and high-order ICA, our analysis suggests that leveraging the TICA can improve our ability to differentiate between this network and others within the default mode domain.

Previous studies show higher model-order ICA results in more spatially granular ICNs (Allen et al., 2011; Iraji et al., 2022); however, the GICA High Order section shows that the DMNs captured by the TICA approach are different from a single-scale ICA with a high model order. Also, more significant group differences were captured using the TICA approach compared with the high-model-order ICA. These results also emphasize the importance of studying functional brain networks in multiple scales using the TICA approach.

In this study, we focus on only two scales with ICA model orders 20 and 4 to capture networks. Further studies can develop our approach in incremental scales to effectively capture ICNs in multiple scales. Furthermore, future studies can explore differences across the different back-reconstruction approaches (Erhardt et al., 2011). Here, we also focused on the voxelwise group comparison of DMN spatial maps. Future studies should focus on the TICA approach using static functional network connectivity and dynamic functional network connectivity on DMNs as well as other brain networks (Iraji et al., 2022).

Brain functional connectivity can occur at different spatial scales, which has been underappreciated. Here, we proposed a TICA approach as a new tool to construct spatial functional hierarchies of the brain and estimate functional sources across multiple spatial scales. Our results on DMN show that our approach can capture replicable comprehensive spatial maps. Moreover, the TICA can detect significant group differences between HC and SZ that are missed using studying brain networks in a single scale, highlighting that our proposed approach represents a promising new tool for studying functional sources.

Data collection was supported by the National Center for Research Resources at the National Institutes of Health (grants NIH 1 U24 RR021992 and NIH 1 U24 RR025736-01).

Shiva Mirzaeian: Conceptualization; Methodology; Visualization; Writing – original draft. Ashkan Faghiri: Writing – review & editing. Vince D. Calhoun: Conceptualization; Investigation; Supervision; Writing – review & editing. Armin Iraji: Conceptualization; Investigation; Supervision; Writing – review & editing.

Vince D. Calhoun, National Institute of Mental Health (https://dx.doi.org/10.13039/100000025), Award ID: R01MH123610. Vince D. Calhoun, National Science Foundation (https://dx.doi.org/10.13039/100000001), Award ID: 2112455.

Intrinsic connectivity network (ICN):

An estimation of a functional source obtained using ICA and typically referring to the components deemed non-artifactual.

Functional source:

A spatial locality with a temporally synchronized pattern. All voxels within a functional source have similar temporal patterns.

Independent component analysis (ICA):

A multivariate technique that decomposes data into a set of maximally independent components.