Leading basic modes of spontaneous activity drive individual functional connectivity organization in the resting human brain – Communications Biology

Participants and study design

We employed three datasets of R-fMRI data from healthy young adults. The first dataset consisted of multiband R-fMRI data from 970 participants from the publicly available S900 data release of the Human Connectome Project (HCP)⁴⁹. The scanning protocol of the HCP dataset was approved by the Institutional Review Board at Washington University. These subjects underwent repeated R-fMRI runs in two sessions. The second dataset, named the sleep deprivation dataset, included repeated R-fMRI data from 20 participants scanned separately during rested wakefulness and after sleep deprivation⁵⁰. The research was approved by the Institutional Review Board of the Institute of Biophysics (Chinese Academy of Sciences). The third dataset, named the Beijing Zang dataset, included R-fMRI data from 198 participants selected from the 1000 Functional Connectomes Project⁵¹. The research was approved by the Institutional Review Board of the State Key Laboratory of Cognitive Neuroscience and Learning at Beijing Normal University. Informed consent was obtained from all participants of the three datasets. All ethical regulations relevant to human research participants were followed. The first two datasets were used for the main analysis, and the third dataset was used for the replication analysis.

Data acquisition

In the HCP dataset, all participants underwent multimodal MRI scanning with a customized 32-channel SIEMENS 3 T Connectome Skyra scanner at Washington University, USA. Four multiband R-fMRI runs were acquired in two sessions for each participant. Briefly, each session consisted of two runs that were separate phases encoded in the left-to-right and right-to-left directions. The R-fMRI scans were obtained using a multiband gradient-echo-planar imaging sequence (repetition time [TR] = 720 ms and 1200 volumes per run, i.e., 14.4 min), with participants’ eyes fixated on a bright projected crosshair. Here, we used only the left-to-right-encoded scans to reduce the potential influence of the phase-encoding directions^15,85. In the original S900 data release, 837 participants completed the left-to-right-encoded R-fMRI scans in both sessions, denoted as REST1 and REST2 separately. Of these, 137 participants were excluded due to missing time points (N = 10), excessive head motion (N = 105) (see “Data Preprocessing”), and arachnoid cysts (N = 22). Data from the remaining 700 participants (aged 21–35 years, M/F: 304/396) were used for the main analyses.

In the sleep deprivation dataset, 20 participants underwent repeated R-fMRI scans separately during rested wakefulness and after sleep deprivation (for design details, see Zhou et al.⁵⁰). All the participants were right-handed and had no history of neuropsychiatric disorders. The MRI data were acquired using a 64-channel 3 T Siemens Prisma scanner at the Beijing MRI Center for Brain Research of the Chinese Academy of Sciences. R-fMRI data were acquired using a T2^*-weighted gradient-echo-planar imaging sequence (TR = 1000 ms and 480 volumes per run), with participants’ eyes fixated on a crosshair. Structural images were acquired using a 3D T1-weighted, magnetization-prepared rapid acquisition gradient-echo sequence. One participant was excluded due to excessive head motion (see “Data Preprocessing”). Data from the remaining 19 participants (aged 18–26 years, M/F: 7/12) were used for the main analysis.

For the Beijing Zang dataset, 198 participants underwent MRI scanning using a 12-channel Siemens Trio Tim 3.0 T scanner in the Imaging Center for Brain Research, Beijing Normal University. R-fMRI data were acquired with participants’ eyes closed (TR = 2000 ms and 235 volumes). One participant was excluded due to differences in scanning orientation, leaving 197 participants (aged 18–26 years, M/F: 75/122) used for the cross-validation analysis.

Data preprocessing

For the HCP dataset, we employed the minimally preprocessed R-fMRI data⁸⁶, followed by ICA-Fix denoising⁸⁷. Four additional steps were performed using the GRETNA package⁸⁸, including the removal of the first 10-second volumes (i.e., 15 volumes), linear detrending, nuisance regression, and temporal filtering (0.01–0.08 Hz). During the nuisance regression, white matter, cerebrospinal fluid, and global signals were included as regressors to further remove the influence of head motion and physiological noise⁸⁹.

The sleep-deprivation dataset and the Beijing Zang dataset were preprocessed with the same pipeline using the GRETNA package⁸⁸. Specifically, the preprocessing included the removal of the first 10-s volumes, realignment, spatial normalization to the Montreal Neurological Institute (MNI) space with the T1-unified segmentation algorithm⁹⁰, linear detrending, nuisance regression, and temporal filtering (0.01–0.08 Hz). During the nuisance regression, we included Friston’s 24 head-motion parameters⁹¹, white matter, cerebrospinal fluid, and global signals as regressors to reduce the influence of head motion and physiological noise⁸⁹.

For these three datasets, we excluded participants with excessive head motion in any scan, including a translation/rotation greater than 3 mm or 3° and a mean framewise displacement (FD) over time⁹² greater than 0.5 mm. After applying these criteria, 105 participants were excluded from the HCP dataset, and 1 participant was excluded from the sleep-deprivation data.

Eigen-microstate analysis of spontaneous brain activity

We applied a recently proposed eigen-microstate analysis^47,48 from the statistical physics theory to the HCP dataset to identify basic modes underlying the rich repertoire of brain activity. The eigen-microstate analysis is an approach specifically designed to extract basic activity modes (i.e., eigen-microstates) of the temporal evolution of a system. This approach has been applied to several complex systems (e.g., the Earth system and stock markets) and reveals meaningful activity modes that show well-defined spatial patterns⁴⁷. In this study, the eigen-microstate analysis was applied to the time courses of each R-fMRI run (i.e., REST1 and REST2) separately. The brain microstates were defined at the nodal level to reduce computational burden.

Definition of the ensemble matrix and microstates

First, we defined 1000 cortical nodes based on a prior functional parcellation⁵², which would enhance functional homogeneity within each nodal region. Then, we extracted the time courses of these nodes for each participant. The time course for each node was further transformed into z-score values with zero mean and unit variance over time. Finally, the normalized nodal time courses were concatenated across participants, resulting in an N × M time course matrix A, where N denotes the number of nodes (i.e., 1000 here) and M denotes the number of time points in the concatenated time courses (i.e., 1185 × 700). Matrix A was considered as an ensemble matrix representing a rich repertoire of brain activity, and each column, A_t, represents a microstate of brain activity at a specific time point from a statistical physics perspective.

Theoretical derivation of eigen-microstates

The eigen-microstate analysis assumes that microstates may exhibit commonalities over time, with instantaneous brain activity arising from the combination of multiple basic modes. It aims to identify the basic modes (i.e., eigen-microstates) that represent fundamental building blocks of rich microstates and are independent of each other^47,48. Specifically, we first calculated the covariance matrix between microstates as C = A^TA⁴⁷, which reflects the spatial similarity between all microstates. Next, we computed the eigenvectors of the matrix C,

wherein v_i is the ith eigenvector of the matrix C, and λ_i is the corresponding eigenvalue. The representativeness of spatial patterns of the microstates at different time points is captured by the eigenvectors⁹³, which take into account the entire spatial similarity pattern between all microstates. Finally, we obtained each eigen-microstate as the weighted sum of the origin microstates:^47,48

wherein v_ti denotes the tth element of eigenvector v_i and quantifies the contribution of tth microstate to the ith eigenvector. In other words, the basic modes (i.e., eigen-microstates) were identified by incorporating brain activity patterns over time and thus can be intuitively interpreted as typical activity patterns. The weight factor of a given eigen-microstate E_i is defined as w_i = E_i^TE_i = λ_i.

Calculation of eigen-microstates. In practical analysis, it is difficult to derive the eigenvectors of the covariance matrix C one by one due to the large number of time points M. Therefore, we obtained the eigen-microstates with the aid of the singular value decomposition (SVD)⁴⁷. Based on the SVD, the ensemble matrix A_N × M was factorized as the product of three matrices:

wherein U = [u₁, u₂, …, u_N] and V = [v₁, v₂, …, v_M] contain eigenvectors for two matrices of AA^T and A^TA, respectively, and Σ_N×M is a diagonal matrix of singular values (σ_i). By substituting this decomposition into Eq. (2), the definition of the eigen-microstate can be reformulated as:

wherein u_i is the eigenvector of the matrix AA^T. The weight of each eigen-microstate can be obtained as w_i = E_i^TE_i = σ_i². The matrix A was normalized by dividing the root-sum-square of all its elements (i.e., a dataset-dependent constant S) before SVD to ensure Σσ_i² = 1.

Of note, the estimation for eigen-microstates in Eq. (4) is similar to the PCA⁹⁴. However, the eigen-microstate analysis and PCA are used in different contexts. The former is specific to the temporal evolution of a system and is used to identify the typical spatial patterns of microstates by analyzing the spatial relationship between time points (see Eq. (2)). The latter, however, is a technique commonly used to reduce the dimensionality of a dataset by finding the principal components that capture the main variance in the multivariate data. In this sense, PCA focuses on the relationship between data variables (here, brain regions) rather than the relationship between time points. Moreover, compared to PCA, which only captures the relative weights of variables, eigen-microstate analysis can provide an intuitive biological interpretation of the eigen-microstates as typical activity states (see Eq. (2)).

To determine whether the brain is dominated by a small number of basic modes (i.e., low-dimensional representations), we identified leading basic modes, whose weights should be⁹⁵: (i) substantially greater than the weight of the subsequent basic mode, known as Cattell’s scree test⁹⁶; (ii) greater than the average weight across all N possible basic modes, i.e., 1/N; (iii) statistically significant according to a permutation test. In each permutation instance, the labels of nodal regions at each time point were randomly shuffled to disrupt the spatial organization. The statistical significance level of each leading basic mode was determined by comparing its original weight (i.e., σ_i²) to the null distribution of the corresponding weights obtained from the 10,000 permutation instances. In our analysis, Cattell’s scree test was performed by identifying the elbow point on the weight curve according to the Kneedle algorithm⁹⁷ (https://github.com/arvkevi/kneed).

We further investigated the spatial patterns of the leading basic modes based on a prior functional system definition⁵³. Seven systems were considered, including the default-mode network, frontoparietal network, limbic network, ventral attention network, dorsal attention network, somatomotor network, and visual network. For each leading basic mode, we estimated the mean fluctuation amplitude for each system by averaging the nodal values within this system.

We also performed the eigen-microstate analysis for each participant to investigate the presence of the leading basic modes at the individual level. In this condition, matrix A in the above analysis was replaced as the normalized time course within each participant for each R-fMRI run.

Cognitive function associations of the leading basic modes

We investigated the potential functional roles of the leading basic modes from two perspectives. First, we examined the association between these leading basic modes and cognitive functions based on the NeuroSynth meta-analytic database (www.neurosynth.org)⁵⁵. For each leading basic mode, we calculated its spatial similarity with all available meta-analytic activation maps using Pearson’s correlation across voxels. The associated cognitive terms are illustrated using word cloud plots.

Second, we compared each of the leading basic modes with 12 cognitive components⁵⁶. Each cognitive component represents a basic activation probability map that is involved in various cognitive tasks⁵⁶. For each cognitive component, we estimated the corresponding node-level version by averaging the activation probabilities of all voxels within each node. We then calculated the spatial similarity between each of the leading basic modes and the 12 cognitive components by using Pearson’s correlation across nodes. To correct for the potential influence of spatial autocorrelation, the statistical significance of each spatial similarity was tested using a permutation test (n = 10,000). The significance level was determined by comparing the original similarity to the null distribution of the corresponding similarity obtained from 10,000 permutation instances. For each permutation instance, we generated a surrogate basic mode map that preserved the spatial autocorrelation of the original basic mode⁵⁷.

Relationship between leading basic modes and FC

Since the leading basic modes dominated the spontaneous fluctuations of brain activity, we further investigated how they contribute to the FC between regions.

The original FC between two nodal regions is defined as Pearson’s correlation between their time courses⁴. As each regional time course (A_it, t = 1, … M) was normalized over time (i.e., mean = 0 and SD = 1), the FC between nodes i and j can be estimated as:^13,43

where M denotes the number of time points in the time course.

By substituting Eq. (3) and Eq. (4) into Eq. (5) and considering the time independence between basic modes, we found FC_ij between nodes i and j can be rewritten as:

where N is the number of all possible basic modes, and E_ik is the ith element of the kth basic mode. Thus, the FC between the two nodes can be attributed to the joint contribution of their coactivation patterns in each basic mode. Notably, the FC_ij estimated from Eq. (6) should be further multiplied by a constant S² to correct for the normalization effect of matrix A prior to the SVD analysis.

To validate the effectiveness of the above theoretical model (i.e., Eq. (6)), we reconstructed the FC matrix according to Eq. (6) by gradually increasing the number of basic modes of interest. We then compared the spatial similarity between the reconstructed and original FC matrices. The spatial similarity was quantified with Pearson’s correlation across the lower triangular elements in the matrices. Specifically, we reconstructed the FC matrix at both the population and individual levels. At the population level, the leading basic modes were obtained from the concatenated normalized time course across all participants. At the individual level, the leading basic modes were identified from the time courses of each participant. We then calculated the similarity between the reconstructed and the original FC matrices for each individual.

We further explored whether the basic modes, especially these leading basic modes, could capture the individual functional organization. First, we estimated the reliability of the reconstructed FC matrix between two runs at the individual level. Given a participant of interest, we evaluated the intra-subject similarity of the reconstructed FC matrices between two runs. We also estimated the inter-subject similarity of this subject as the averaged spatial similarity of this participant in the first run (i.e., REST1) with all the other participants in the second run (i.e., REST2). Next, we examined the individual uniqueness in the reconstructed FC matrices by performing an FC-based individual identification analysis between two runs (i.e., REST1 and REST2)¹³. For each participant, we compared the reconstructed FC matrix of this participant in REST1 with those of all the participants in REST2. If the participant with the highest similarity in REST2 was the same participant given in REST1, the identification was correct; otherwise, it was incorrect. Identification accuracy was defined as the proportion of participants that were correctly identified. A higher accuracy indicates a higher capability of the FC matrix in distinguishing individuals. We also evaluated the differential identifiability⁵⁸ to quantify the individual uniqueness of the reconstructed FC matrices. The differential identifiability, I_diff, was defined as I_diff = (I_self − I_others) × 100, wherein I_self and I_others denoted the mean intra-subject similarity and mean inter-subject similarity between two runs, respectively. A larger I_diff value indicates more individual-specific information is captured. For comparison, individual identification analysis was also performed based on the original FC matrix.

Influence of sleep deprivation on the leading basic modes

To assess whether the leading basic modes are affected by mental states, we applied the eigen-microstate analysis to the sleep deprivation dataset⁵⁰. We identified the leading basic modes (see “Eigen-microstate analysis”) separately from the R-fMRI data obtained in the two states represented in the dataset (i.e., rested wakefulness vs. post-sleep deprivation). First, we examined the spatial correspondence of the leading basic modes determined at rested wakefulness with those obtained from REST1 of the HCP dataset to investigate the reproducibility of the leading basic modes. Next, we compared the basic modes obtained at the different states (i.e., rested wakefulness vs. post-sleep deprivation) to examine the potential influence of sleep deprivation.

For each leading basic mode that maintained spatial correspondence between two states, we tested differences in regional fluctuation amplitudes between the two states by using a permutation test (n = 10,000). In each permutation instance, the state labels of the R-fMRI data were shuffled for each participant. Multiple comparisons across nodal regions were corrected using the FDR approach (corrected p < 0.05)⁹⁸. Given that the basic mode showed significant changes, we further investigated how interregional coactivation patterns differed between the two states. Briefly, we obtained the system-level coactivation pattern for the rested wakefulness and post-sleep deprivation separately. The within-system and between-system coactivation values were obtained by averaging the interregional coactivation values within the same system and between different systems, respectively. Significance levels of differences in the coactivation pattern were also estimated using a permutation test (n = 10,000) and corrected for multiple comparisons (FDR corrected p < 0.05).

Validation analysis

The reliability of the leading basic modes was validated by considering five analysis strategies that may affect the identification of the leading basic modes. In each case, the leading basic modes were re-estimated and compared with those obtained in the main analyses (i.e., REST1 in HCP). (i) Number of participants. We re-identified the leading basic modes by including R-fMRI data from different numbers of participants from the HCP dataset, with the number of participants varying from 50 to 650 with a step of 50. Given a participant number, we randomly sampled participants from all the 700 participants, and this process was repeated 100 times. For each sampling instance, we counted the number of the leading basic modes and compared the spatial patterns of the first five basic modes with those obtained from the whole population. (ii) Head motion. Head motion during R-fMRI scanning can affect the fluctuation amplitudes of BOLD signals⁹⁹. Different from the main analysis, we used stricter head motion exclusion criteria for R-fMRI data in the HCP dataset (i.e., >2 mm or 2° in any direction or mean FD > 0.2 mm) to further reduce the influence of head motion. (iii) Global signal regression. In the main analysis, the global signal was regressed to better reduce the influence of head motion and non-neural signals^89,99. To assess the potential influence of the global signal, we re-preprocessed the R-fMRI data in the HCP dataset without global signal regression. (iv) Brain parcellation. To assess the influence of spatial resolution, we extracted regional time courses from the HCP dataset by using the same type of functional parcellations with different spatial resolutions (i.e., comprising 200 and 400 cortical regions)⁵². The leading basic modes obtained from different spatial resolutions were compared at the functional system level⁵³ and the voxel-wise level. In the latter case, the voxels within the same nodal regions were assigned the same amplitude values for each basic mode, regardless of the spatial resolution. In cases (i)–(iv), the validation analysis was performed based on REST1 of the HCP dataset. (v) Reproducibility across datasets. We identified the leading basic modes from another independent dataset, i.e., the Beijing Zang dataset⁵¹, and compared them with those in the HCP dataset.

Statistics and reproducibility

Permutation tests were performed to obtain the statistical significance of the following analyses, including the number of the leading basic modes (Fig. 1a), spatial similarities between the leading basic modes and cognitive components (Fig. 2b), and the influence of sleep deprivation on the fluctuation amplitudes and coactivation patterns of the leading basic modes (Fig. 4c, d). For each analysis, 10,000 permutation instances were generated during the permutation test. Paired t-tests were performed to test the statistical differences between intra- and inter-subject similarity of the reconstructed FC matrices between two runs (n = 700) (Fig. 3d). Five analysis strategies were considered to verify the reproducibility, including (i) varying the number of participants (n ranged from 50 to 650 with a step of 50); (ii) using stricter head motion exclusion criteria (n = 415); (iii) performing nuisance regression without global signal regression (n = 700); (iv) defining brain nodes based on two functional parcellations with different spatial resolutions (n = 700); and (v) using another independent dataset, i.e., the Beijing Zang dataset (n = 197). The number and spatial patterns of the leading basis modes were examined in these cases.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.