Abstract
The human brain exhibits a high degree of individual variability in both its structure and function, which underlies intersubject differences in cognition and behavior. It was previously shown that functional connectivity is more variable in the heteromodal association cortex but less variable in the unimodal cortices. Structural connectivity (SC) is the anatomical substrate of functional connectivity, but the spatial and temporal patterns of individual variability in SC (IVSC) remain largely unknown. In the present study, we discovered a detailed and robust chart of IVSC obtained by applying diffusion MRI and tractography techniques to 1,724 adults (770 males and 954 females) from multiple imaging datasets. Our results showed that the SC exhibited the highest and lowest variability in the limbic regions and the unimodal sensorimotor regions, respectively. With increased age, higher IVSC was observed across most brain regions. Moreover, the specific spatial distribution of IVSC is related to the cortical laminar differentiation and myelination content. Finally, we proposed a modified ridge regression model to predict individual cognition and generated idiographic brain mapping, which was significantly correlated with the spatial pattern of IVSC. Overall, our findings further contribute to the understanding of the mechanisms of individual variability in brain SC and link to the prediction of individual cognitive function in adult subjects.
Significance Statement
White matter (WM) connectivity between gray matter regions plays an important role in integrating information from distributed regions when individuals perform complex cognitive functions. Unique WM connectivity is the neuroanatomical identity of each individual. The authors systematically explored the spatiotemporal pattern of individual variability in structural connectivity (SC), and the results show higher variability in SC relates to higher neuroplasticity. In addition, this study reveals that SC involved in executive function and attention task is different among individuals and highlights the importance of individualized statistical methods for mapping neural pathway of complex cognition.
Introduction
Much like fingerprints, everyone's brain is unique. The human brain exhibits great individual variability in anatomy, structure, and function (Rademacher, 2001; Bürgel et al., 2006; Sugiura et al., 2007), which likely underlies the individual differences in behavior and cognition. The individual differences among human brains are jointly shaped by genetic and environmental factors (Glendenning and Masterton, 1998). Unraveling the spatial and temporal patterns of individual differences can facilitate our understanding of the neural underpinnings of interindividual differences in behavior and cognition.
The degree of variation in the human brain among different individuals is not uniform across the cortex. For example, functional MRI (fMRI) studies have revealed that individual differences in functional connectivity are significantly higher in the heteromodal association cortex and lower in unimodal cortices. This distinction is significantly correlated with the degree of evolutionary cortical expansion, suggesting a potential evolutionary root of functional variability (Mueller et al., 2013). In terms of cortical morphology, individual variability in cortical folding is higher in association cortex areas than in the motor cortex (Hill et al., 2010). For white matter (WM) tracts, long association fibers, e.g., the superior longitudinal fascicle, are more variable in the position and extent of their fiber tracts than the corticospinal tract and optic tract (Bürgel et al., 2006). These findings consistently demonstrate spatial distributions of individual variability in brain structure and function, which is potentially modulated by the underlying biological, genetic mechanisms and evolutionary pressure.
Advanced brain imaging and network modeling techniques have made it possible to form the concept of the human brain connectome, providing a new perspective to investigate the interaction and integration among different brain regions. Individual variability in functional connectivity (IVFC) has been comprehensively examined with fMRI techniques (Mueller et al., 2013; Gao et al., 2014; Xu et al., 2019; Ma et al., 2021). The IVFC can be used to predict interindividual variation in behavior. It can also be used as a “fingerprint” for individual identification (Finn et al., 2015). Moreover, the IVFC shows dynamic changes with lifespan: it initially emerges in the third trimester, decreases with brain maturation, and increases with the aging process, especially in the heteromodal association cortices. With diffusion MRI (dMRI) and tractography algorithms, WM structural connectivity (SC) among brain regions can be reconstructed in vivo. Unlike IVFC, the spatial and temporal patterns of individual variability in SC (IVSC) remain largely unknown. To date, only three studies have reported spatial distributions of IVSC and showed inconsistent results (Chamberland et al., 2017; Mansour et al., 2021; Sun et al., 2022). The conflicting findings may result from the heterogeneity of the samples, acquisition parameters, and methodological differences. Therefore, a reliable and robust characterization of the spatial and temporal patterns of IVSC is warranted.
Importantly, individual differences in brain structural and functional connectivity can provide sensitive features for individual prediction. Recently, brain imaging features combined with machine learning approaches have been proposed for cognitive prediction at the individual level (Sui et al., 2020). However, traditional machine learning models applying the same features and features’ weights to different individuals are still a population-level predictive approach that relies heavily on multivariate pattern information conserved across individuals. Given the substantial individual differences in brain structure and function (D. Wang et al., 2015; Cui et al., 2020), individualized machine learning methods are particularly important for understanding the neural underpinnings of cognition, as they can capture the complexity of the neural mechanisms and representations of cognition at the individual level.
To address this gap, we first aimed to characterize a reliable and replicable spatiotemporal pattern of IVSC by utilizing multiple imaging datasets. Second, we sought to explore the biological and genetic mechanisms associated with IVSC by integrating several brain atlases into the analysis. Third, we aimed to design a machine learning framework capable of generating idiographic brain mapping to predict cognition and assess the association between IVSC and the variation in idiographic brain mapping.
Materials and Methods
Participants
The human MRI data employed in this work came from five datasets, HCP-YA (Human Connectome Project-Young Adults; Van Essen et al., 2013), HCP Retest (Van Essen et al., 2013), HCP-A (HCP-Aging; Harms et al., 2018), Cam-CAN (Cambridge Centre for Aging and Neuroscience; Cam-CAN et al., 2014), and BABRI (Beijing Aging Brain Rejuvenation Initiative; Yang et al., 2021). All studies were approved by the research ethics review boards of the respective institutions. This project was approved by the research ethics review board of Beijing Normal University. These five datasets are described below.
HCP-YA
We used the final release of HCP-YA, which includes 1,206 healthy young adult participants. We restricted our analysis to subjects from different families to avoid the effects of shared genes and the environment on the evaluation of individual variability (n = 457 subjects remained following this exclusion). After 38 subjects were excluded due to the low quality of the raw diffusion-weighted MRI scans or the data in the QSIPrep pipeline (Cieslak et al., 2021), the present study included data from a total of 419 healthy young adult participants from the HCP-YA study (221 males, age 22–37 years).
HCP Retest
Among the 1,206 subjects from the HCP-YA study, 46 subjects underwent two separate scans with an interval ranging from 0.5 to 11 months. Three subjects were excluded after quality control of MRI scans and the QSIPrep pipeline (Cieslak et al., 2021). The two MRI scans of the remaining 43 subjects from HCP Retest were included in this study (13 males, age 22–35 years).
HCP-A
HCP-A is an ongoing study aimed at enrolling 1,200+ healthy adults (age 36–100+ years). We used HCP-A Release 2.0, which included 725 healthy adults. Considering the influence of gray matter atrophy due to aging on the evaluation of IVSC, we excluded participants who were older than 68 years, which is an upper tertile of ages of all participants from the HCP-A, Cam-CAN, and BABRI. After excluding the subjects whose data quality in the raw MRI data or QSIPrep pipeline (Cieslak et al., 2021) was low, data from 470 healthy adults from the HCP-A study were included in this study (300 males, age 36–68 years).
Cam-CAN
Cam-CAN enrolled 2,681 healthy participants, and 700 of those participants underwent MRI scans. Similarly, participants older than 68 years were excluded. We also excluded participants younger than 30 years because we wanted to focus only on IVSC in the aging process. Finally, 382 participants from Cam-CAN were included in this study (195 males, age 31–68 years).
BABRI
BABRI is an ongoing study based on a community cohort focusing on asymptomatic stages of dementia and aiming to develop prevention strategies for cognitive impairment. Similarly, after excluding participants older than 68 years, 414 healthy participants from the BABRI were included (141 males, age 45–68 years).
Cognitive assessment
To explore the association between cognition and individual variability, we only included the cognitive data from the HCP-A and BABRI because both studies used the same neuropsychological tests to assess participants’ executive function and attention. For the HCP-A dataset, in addition to the NIH toolbox, additional cognitive tests, including the Rey Auditory Verbal Learning Task, delay discounting, and Trials A and B, were administered to all participants as cognitive assessments. For the BABRI dataset, a battery of neuropsychological tests was adopted to assess cognitive function across five domains: memory, executive function, attention, visuospatial ability, and language. Among the five domains, executive function was measured via the Trail B, Stroop B and C tests; attention was assessed by the Trail A, symbol digit modalities and Stroop A tests; and memory was examined by the Rey Auditory Verbal Learning Task.
Image acquisition
HCP-YA and HCP Retest
MRI data were collected with a customized 3 T Connectome Scanner adapted from Siemens Healthineers Skyra. T1-weighted (T1w) scans used 3D magnetization-prepared rapid gradient echo (MPRAGE; slices, 256; TR, 2,400 ms; TE, 2.14 ms; flip angle, 8°; voxel size, 0.7 mm isotropic). A multishell diffusion-weighted echoplanar imaging sequence was used for diffusion-weighted MRI (dMRI) data (90 diffusion-weighted directions for b = 1,000, 2,000 and 3,000 s/mm2 and 18 images with b = 0 s/mm2; slices, 111; TR, 5,500 ms; TE, 89.50 ms; slice thickness, 1.25 mm; FOV, 210 × 180 mm2; and acquisition matrix, 168 × 144).
HCP-A
MRI data were acquired with a 3 T MRI scanner (Prisma, Siemens Healthineers). T1w scans used multiecho MPRAGE (sagittal slices, 208; TR, 2,500 ms; TE, 1.81/3.6/5.39/7.18 ms; TI, 1,000 ms; slice thickness, 0.8 mm; flip angle, 8°; FOV, 256 × 256 mm2; and acquisition matrix, 320 × 320). A multishell diffusion-weighted echoplanar imaging sequence was used for dMRI data (92 diffusion-weighted directions for b = 1,500; 93 diffusion-weighted directions for 3,000 s/mm2; 14 images with b = 0 s/mm2; slices, 92; TR, 3,230 ms; TE, 89.20 ms; slice thickness, 1.5 mm; FOV, 210 × 210 mm2; and acquisition matrix, 140 × 140).
Cam-CAN
All MRI datasets were collected at a single site using a 3 T MRI scanner (Trio, Siemens Healthineers). 3D MPRAGE was used for the T1w scan (TR, 2,250 ms; TE, 2.99 ms; TI, 900 ms; flip angle, 9°; FOV, 256 × 240 × 192 mm; voxel size, 1 mm isotropic). A twice-refocused spin‒echo sequence was used for dMRI (30 diffusion gradient directions for b = 1,000 and 2,000 s/mm2; three images with b = 0 s/mm2; axial slices, 66; TR, 9,100 ms; TE, 104 ms; FOV, 192 × 192 mm).
BABRI
MRI data were acquired with a 3 T MRI scanner (Trio, Siemens Healthineers MAGNETOM). A sagittal 3D MPRAGE was used for the T1w scan (sagittal slices, 176; TR, 1,900 ms; TE, 3.44 ms; TI, 900 ms; slice thickness, 1 mm; flip angle, 9°; FOV, 256 × 256 mm2; and acquisition matrix, 256 × 256). A single-shell diffusion–weighted echoplanar imaging sequence was used for dMRI data (30 diffusion-weighted directions with b = 1,000 s/mm2 and a single image with b = 0 s/mm2; axial slices, 70; TR, 9,500 ms; TE, 92 ms; slice thickness, 2 mm; FOV, 256 × 256 mm2; and acquisition matrix, 128 × 128).
Image processing
The image processing pipeline consisted of preprocessing and reconstruction workflows. QSIPrep (Cieslak et al., 2021) 0.14.3, which is based on FSL (Andersson and Sotiropoulos, 2016), DSI studio (Yeh and Tseng, 2011), DIPY (Garyfallidis et al., 2014), ANTs (Avants et al., 2008), and MRtrix (Tournier et al., 2019), was used to perform image processing and create individual WM networks. Detailed procedures are described below.
Preprocessing
For T1w images, ANTs was used to extract the brain data and generate a brain mask. FSL was utilized for brain tissue segmentation. Finally, ANTs was used to normalize the T1w image to the MNI space. For dMRI data, images were grouped according to the phase encoding direction. Then, MRtrix was used to denoise diffusion-weighted images. FSL was used for distortion and head motion correction. Finally, the template image with b = 0 was created and registered to the skull-stripped T1w image.
Reconstruction
The reconstruction workflow performed advanced reconstruction and tractography methods to create a WM network. Single-shell and multishell-constrained spherical deconvolutions were separately applied to single-shell and multishell dMRI data to estimate voxel-level diffusion models. A probabilistic tractography algorithm in MRtrix was used to produce a fixed number of streamlines, and an anatomically constrained tractography framework was used during tracking. Then, spherical deconvolution-informed filtering of tractograms was applied to improve the quantitative nature of whole-brain streamline reconstructions. The Brainnetome Atlas was generated using a connectivity-based parcellation framework, demonstrating its ability to represent human connectome (Fan et al., 2016). There, it was utilized to define nodes in the WM network, with the number of streamlines connecting a pair of regions used to weight edges.
Analyses of IVSC
Harmonization of multisite data
For removing the site-specific effects in SC obtained from multiple sites, the ComBat model was utilized. ComBat was previously applied successfully for harmonizing measures from MRI data, such as cortical thickness (Fortin et al., 2018) and functional connectivity (A. A. Chen et al., 2022). This approach assumes that sites have both additive and multiplicative effects on data. Hence, the model is defined as follows:
Definition of individual variability of the structural network
The IVSC was quantified according to the definition in previous studies (Mueller et al., 2013; Xu et al., 2019; Ma et al., 2021). For a given node
i, the intersubject difference was defined as follows:
The three middle-aged datasets do not include retest data, but different studies have identified a similar spatial distribution of SC reproducibility (Bassett et al., 2011; Buchanan et al., 2014; Zhao et al., 2015). This indicates that the spatial distribution of intrasubject differences is relatively stable. Therefore, intrasubject difference was estimated using the HCP Retest dataset, where each participant underwent two MRI scans, and was applied to other datasets to calculate IVSC. It was defined as follows:
Sliding window analysis
To examine the structural network variability as a function of age, cross-subject sliding window analyses were performed for participants from HCP-A, Cam-CAN, and BABRI in ascending order of age, with a window length of 10 participants and a step size of 1 participant. Using this method, 1,253 overlapping subgroups were generated, and their structural network variabilities and mean ages were calculated.
Spatial permutation testing (Moran spectral randomization)
To evaluate the association between the spatial distribution of structural network variability and other cerebral properties, we used a spatial permutation procedure called Moran spectral randomization (Wagner and Dray, 2015) to control the effect of spatial autocorrelation on the association. The inverse Euclidean distance between nodes was used to denote the spatial proximity when using Moran spectral randomization to generate 10,000 random brain maps that preserved the spatial autocorrelation of the original map. The null distribution of Spearman's r values was obtained by relating structural network variability with random brain maps. The p value was computed as the proportion of null r values that were greater than the true r value.
Transcriptional analyses
Transcriptional data preprocessing
To explore the genetic mechanisms underlying the specific distribution of IVSC, we employed transcriptional data from six postmortem brains provided by the Allen Human Brain Atlas (AHBA; Hawrylycz et al., 2012). Analyses were limited to the left brain only because only two of six brains included samples from the right hemisphere. The Abagen toolbox (Markello and Arnatkevi, n.d.) was used to preprocess the transcriptional data according to the suggestions in a previous study (Arnatkeviciute et al., 2022).
Specifically, microarray probes were reannotated to corresponding genes based on the database provided by Arnatkeviciute et al. (2019). We filtered noisy probes that did not exceed the background in at least 50% of all samples across all subjects. For each gene, the representative probe with the expression level most correlated to RNA-seq measures collected in two of six brains was selected. Next, tissue samples were assigned to the closest regions in the Brainnetome Atlas, while samples >2 mm from all regions were not assigned to any regions. To account for interindividual variation, expression values within each brain across regions were normalized with a scaled robust sigmoid function. To retain the genes with consistent expression patterns across six brains, differential stability (Hawrylycz et al., 2015) was calculated for each gene, and a threshold of 0.1 was imposed. Finally, regional expression values were averaged across donors.
It is important to provide gene specificity on the observed transcriptomic-neuroimaging effect to avoid nonspecific effects (Wei et al., 2021). To retain brain-specific genes, we selected genes with higher expression in brain sites than in nonbrain sites using the Genotype-Tissue Expression database (https://www.gtexportal.org/). This procedure produced a matrix with 123 rows corresponding to brain regions and 2,112 columns corresponding to the retained genes.
PLS analysis
The correlation between gene expression and structural network variability was examined with a partial least square (PLS) analysis which is an unsupervised multivariate method to project two matrices into a lower-dimensional space where the covariance between them is maximal. This analysis was performed using the PLSRegression function in the Python package scikit-learn (https://scikit-learn.org/stable/index.html) with the number of components set as 1 because IVSC only has one dimension.
Bootstrap resampling
To obtain a reliable estimation of gene weights, bootstrap resampling was used. The values in the structural network variability vector were randomly selected with replacement 10,000 times. PLS analysis was performed using these bootstrapped structural network variability vectors resulting in distributions of gene weights. The reliable gene weights were quantified by the ratio of the true gene weight to the standard error of its bootstrapped distribution. The ratio is the combination of both the value and reliability of gene weights.
Gene set enrichment analysis
To further investigate biological pathways related to the genes identified by PLS analysis, the method used by Hansen et al. (2021) was adopted. The latest gene annotation provided by Gene Ontology (geneontology.org) was downloaded on July 2, 2022. For each biological process in gene annotation, the enrichment score was defined as the mean weights of the genes in the biological process. The null distribution of enrichment scores was constructed by permuting the structural network variability map 10,000 times with Moran spectral randomization, repeating PLS analyses, and calculating the enrichment score. This analysis was separately performed on the positive genes with 50% of the most positive weights and negative genes with 50% of the most negative weights. For positive genes, significance was measured by the proportion of null enrichment scores that were greater than the true enrichment score. For negative genes, the proportion of null enrichment scores that were smaller than the true enrichment score was utilized to quantify significance.
To determine cell types where the genes identified by PLS are preferentially expressed, the method performed by Hansen et al. (2021) was adopted. Briefly, gene annotation about the cell type was provided by Seidlitz et al. (2020), who used hierarchical clustering based on spatial expression profiles to distinguish seven canonical cell types: astrocytes, endothelial cells, microglia, excitatory neurons, inhibitory neurons, oligodendrocytes, and oligodendrocyte precursors. For each cell type, the enrichment score was defined as the proportion of genes that were identified by PLS analysis and preferentially expressed in that type as well. This analysis was separately performed on the positive genes with 50% of the most positive weights and negative genes with 50% of the most negative weights. A nonparametric permutation test was used to quantify the significance.
Development of gene expression
We used BrainSpan (Miller et al., 2014), a database of gene expression across development, to depict the development of the gene expression component determined by PLS analysis in the AHBA. Postmortem brains in BrainSpan were grouped into five developmental stages: fetus (8–37 postconception weeks), infant (4 months to 1 year), child (2–8 years), adolescent (11–19 years), and adult (21–40 years). For each group, gene expression scores in 16 regions were calculated by multiplying the gene weights derived from the PLS analysis in the AHBA and mean expression across postmortem brains in the group.
To relate gene expression scores in BrainSpan and structural network variability, we mapped the Brainnetome Atlas into 14 Brodmann areas that are included in BrainSpan. The other two Brodmann areas were not included in subsequent analyses because no area in Brainnetome Atlas was mapped to the two Brodmann areas. For each Brodmann area, structural network variability was averaged across the Brainnetome Atlas regions corresponding to the Brodmann area. Then, the relationship between gene expression scores and structural network variability was examined in each developmental stage.
Cognitive prediction
Considering that the variability in SC, the SC involved in a cognitive task may differ across individuals. To test this hypothesis, we modified the traditional ridge regression to include individualized feature selection by integrating the rectified linear unit activation. We then investigated whether incorporating individualized feature selection into the traditional machine learning method could improve the accuracy of cognitive prediction by comparing the performance of traditional and modified ridge regression (Extended Data Fig. 5-1A). The mean absolute error (MAE) and Pearson's r value were used to evaluate the models’ performance. The thresholds of 0.66 and 0.33 are commonly used (Zeighami et al., 2019) for group thresholding. So, we categorized the structural connections into two groups: consistent edges, which were selected by the modified ridge regression to predict the cognition of at least two-thirds of participants, and individualized edges, which were selected by the modified ridge regression to predict the cognition of less than one-third of participants (Extended Data Fig. 5-2B). We then quantified the number of consistent edges and individualized edges for each region and performed Pearson's correlation between the number of these edges and the IVSC. This was down to explore whether the edges of a region with high IVSC are only selected to predict cognition for a few individuals. Notably, we excluded participants with scores <10 or >300 when predicting executive function and those with scores <10 when predicting attention, as these participants’ cognitive assessments may not have been reliable.
Modified ridge regression
To test whether the SC involved in cognitive tasks differed across individuals, we used two models to predict individual cognition. One was ridge regression, which is formalized as follows:
Considering that the SC encoding cognitions may be different across individuals, we proposed a modified ridge regression that is a neural network with a hidden layer. The model can be defined as follows:
Prediction framework
To train models for cognitive prediction, a nested 2F-CV (twofold cross-validation) where the outer 2F-CV was used in the evaluation of models’ performance and the inner 2F-CV was used to determine the optimal hyperparameter α.
Outer 2F-CV
For the outer 2F-CV, we merged the HCP-A and BABRI datasets because they both used Trail B and Trail A to measure individuals’ executive function and attention, respectively. Then, participants in the merged dataset were sorted according to the outcome (i.e., executive function or attention). The participants with an odd rank were assigned to Subset 1, and other participants were assigned to Subset 2. In the first fold, Subset 1 was used as the training dataset, and Subset 2 was used as the test dataset. In the second fold, Subsets 2 and 1 were used as the training and testing datasets, respectively. Each feature was normalized across the training dataset, and the mean value and standard deviation derived from the training dataset were also applied to the test dataset. After the optimal hyperparameter α was selected, a model with the optimal α was trained using data from all participants in the training dataset and then applied to the test dataset to predict the outcome. Pearson's correlation coefficients and MAE between the predicted and true outcomes of participants in the test dataset were calculated to quantify the prediction accuracy.
Inner 2F-CV
For the inner 2F-CV, the training dataset was further split into two subsets as for the outer 2F-CV. These subsets were used to train models under different values of hyperparameter α ([2−4, 2−3, …, 24, 25]) in turn, and these models were applied to the remaining dataset to calculate Pearson's correlation coefficient and MAE between the predicted and true outcomes. For each value of α, Pearson's correlation coefficients and MAE were averaged across two inner loops. Then, the sum of the normalized mean Pearson's correlation coefficient and the normalized reciprocal of the mean MAE was defined as the inner prediction accuracy. The value of hyperparameter α with the highest inner prediction accuracy was selected as the optimal α.
Randomly split 2F-CV
To validate whether the result is influenced by data splitting, we randomly split the merged dataset into two halves in both outer 2F-CV and inner 2F-CV and calculated the mean Pearson's correlation coefficient and MAE across two outer folds. This procedure was repeated 100 times to obtain the distribution of mean Pearson's correlation coefficients and MAE.
Significance of prediction accuracy
A permutation test was performed to evaluate the significance of the prediction accuracy. Briefly, we performed nested 2F-CV 1,000 times with permuting outcomes across the training dataset in each run. Null distributions of Pearson's correlation coefficients and MAE were obtained. The significance of Pearson's correlation coefficient was determined as the proportion of the null coefficient higher than the true coefficient, and the significance of MAE was defined as the proportion of null values lower than the true value.
Replication analyses
Effects of network construction method
IVSC maps were recalculated by changing one setting of network construction method each time. These settings included whether to use spherical deconvolution-informed filtering of tractograms and different weighting methods—either number of streamlines or the density of streamlines (normalized by average volume of a pair of interconnected regions). The Spearman's correlation between the IVSC map in the main results and these recalculated IVSC maps was used to evaluate the robustness of an IVSC pattern.
Effects of edge threshold
To evaluate whether the IVSC pattern identified in this study is robust to edge threshold, we used different edge thresholds to remove weak edges in structural networks and then recalculated the IVSC maps. The Spearman's correlations between the IVSC map without removing edges and the IVSC maps with removing weak edges were used to assess the stability of IVSC pattern.
Effects of brain parcellation
In addition to the Brainnetome Atlas, we also used two other atlases to define network nodes: the first version of Automated Anatomical Labelling (AAL1) atlas, which is coarse and parcellates the cerebrum into 90 regions, and Gordon's atlas, which is fine and parcellates the cerebrum into 333 regions. The IVSC maps were separately calculated from the networks based on AAL1 and Gordon's atlases and analyzed with the procedure similar to that performed in the networks based on the Brainnetome Atlas.
Data availability
HCP-YA, HCP Retest, and HCP-A are available at https://www.humanconnectome.org/. Cam-CAN is available at https://www.cam-can.org/index.php?content=dataset. BABRI is available upon a request. AHBA is available at https://help.brain-map.org/display/humanbrain/Allen+Human+Brain+Atlas. The Genotype-Tissue Expression database is available at https://www.gtexportal.org/.
Code accessibility
The code for generating individual variability in SC, statistics, and cognitive prediction is available on GitHub (https://github.com/forwho/Individual-Difference-in-Structural-connectivity).
Results
Spatiotemporal pattern of IVSC
To depict the spatiotemporal pattern of IVSC, 1,724 participants from the HCP-YA (Van Essen et al., 2013), HCP Retest (Van Essen et al., 2013), HCP-A (Harms et al., 2018), Cam-CAN (Cam-CAN et al., 2014), and BABRI (Yang et al., 2021) were studied. Participants’ demographic information is shown in Table 1. Trail A and Trail B tests were used to evaluate the attention and executive function of the participants, respectively.
Demographical information of participants from the five datasets
The HCP Retest, where participants underwent two MR scans within 1 year, was used to estimate intraindividual differences in SC and regional regression coefficients between interindividual and intraindividual differences in SC. A small intraindividual difference in SC indicated high reliability of SC (Extended Data Fig. 1-1). The intraindividual difference and regression coefficients were applied to the other datasets to estimate IVSC.
To obtain the IVSC map, we merged three middle-aged datasets—HCP-A, Cam-CAN, and BABRI—after removing site effect with ComBat (Fortin et al., 2018). As shown in Figure 1B, the limbic network exhibits the highest variability, association networks show moderate variability, and unimodal networks display the lowest variability. To assess the stability of the IVSC maps in middle-aged adults across different datasets, we calculated IVSC maps of HCA-A, Cam-CAN, and BABRI individually (Extended Data Fig. 1-2A–C). The results indicated that the IVSC maps derived from these datasets were not only significantly correlated to each other but were also with the IVSC map derived from the combined dataset (Fig. 1C; all Spearman's r(246) values >0.74). To evaluate the effect of age on the spatial pattern of IVSC in adults, we also calculated the IVSC of young adults (Fig. 1A) from the HCP-YA. The IVSC map of young adults was similar to that of middle-aged adults from the combined dataset (Spearman's r(246) = 0.81; p < 0.001). In addition, we split the combined dataset into age tertiles and calculated IVSC for each age group. The three IVSC maps were highly similar to each other (Extended Data Fig. 1-3; all Spearman's r(246) values ≥0.97). In addition, by using sliding window analysis, we found that IVSC significantly increased with age in almost all brain regions (Fig. 1D). Both the mean IVSC and standard deviation of IVSC were positively correlated with age (Fig. 1E,F).
Spatiotemporal pattern of IVSC. A, The IVSC across participants in HCP-YA. B, The IVSC across participants in the combination of HCP-A, Cam-CAN, and BABRI studies. Both IVSCs were corrected for the intrasubject difference shown in Extended Data Figure 1-1. C, The correlation between IVSC values from different datasets. D, Age-related rates of IVSC obtained by performing sliding window analysis in the combined dataset of HCP-A, Cam-CAN, and BABRI. E, The relationship between mean IVSC and mean age. F, The relationship between the standard deviation of IVSC and mean age. We evaluated the effect of age on the IVSC and the impact of different processing pipelines on the IVSC. The IVSCs derived from each middle-aged dataset were shown in Extended Data Figure 1-2. The IVSCs in three age-based tertiles were shown in Extended Data Figure 1-3. The IVSCs based on different network construction methods were shown in Extended Data Figure 1-4. The IVSCs under different edge thresholds were shown in Extended Data Figure 1-5. The IVSCs using different atlases were shown in Extended Data Figure 1-6. IVSC, individual variability in structural connectivity; HCP-YA, Human Connectome Project-Young Adults; HCP-A, Human Connectome Project Aging, Cam-CAN, Cambridge Centre for Aging and Neuroscience; BABRI, Beijing Aging Brain Rejuvenation Initiative; SD, standard deviation.
Figure 1-1
Intraindividual differences in structural connectivity. The intraindividual differences in structural connectivity were calculated based on HCP Rest datasets. HCP, Human Connectome Project. Download Figure 1-1, TIF file.
Figure 1-2
Individual variability in structural connectivity maps derived from three different datasets. A, Individual variability in structural connectivity maps derived from the HCP-A. B, Individual variability in structural connectivity maps derived from the Cam-CAN. C, Individual variability in structural connectivity maps derived from the BABRI. HCP-A, Human Connectome Project Ageing, Cam-CAN, Cambridge Centre for Ageing and Neuroscience; BABRI, Beijing Ageing Brain Rejuvenation Initiative. Download Figure 1-2, TIF file.
Figure 1-3
Individual variability in structural connectivity in three age-based tertiles. We split the samples in the combined dataset into age tertiles and calculated individual variability in structural connectivity in each group. We found that individual variability in structural connectivity was nearly identical in each of the three age-based tertiles. Download Figure 1-3, TIF file.
Figure 1-4
Individual variability in structural connectivity based on different network construction methods. A, Individual variability in structural connectivity without spherical deconvolution-informed filtering of tractograms. B, Individual variability in structural connectivity based on the network weighted by the density of streamlines. Download Figure 1-4, TIF file.
Figure 1-5
Effects of edge thresholds on IVSC measurement. We used different edge thresholds to remove weak edges in structural networks and then re-calculated the IVSC maps. We found that the IVSC maps after removing weak edges were highly similar to the IVSC map without removing edges. Inset plots show the IVSC maps with removing weak edges under different threshold. IVSC, individual variability in structural connectivity. Download Figure 1-5, TIF file.
Figure 1-6
Temporal spatial pattern of IVSC based on AAL and Gordon’s atlases. A, The IVSC map based on AAL atlas. B, Age-related rates of IVSC based on AAL atlas. C, the relationship between mean IVSC based on AAL atlas and mean age. D, The relationship between the standard deviation of IVSC based on AAL atlas and mean age. E, The IVSC map based on Gordon’s atlas. F, Age-related rates of IVSC based on Gordon’s atlas. G, the relationship between mean IVSC based on Gordon’s atlas and mean age. H, The relationship between the standard deviation of IVSC based on Gordon’s atlas and mean age. IVSC, individual variability in structural connectivity; AAL, automated anatomical labelling. Download Figure 1-6, TIF file.
IVSC was also assessed within eight brain systems, including seven cortical networks (Yeo et al., 2011; Fig. 2A) and a subcortical system. The limbic network showed the highest level of IVSC, whereas the somatomotor network and visual network were the least variable (Fig. 2B). A sliding window analysis revealed that the mean variability in all eight brain systems significantly increased with age (Fig. 2C–J).
IVSC quantified in eight brain systems. A, The spatial distribution of seven cortical networks defined by Yeo and his colleagues (Yeo et al., 2011). B, Mean IVSC in eight brain systems. C–J, The relationship between mean IVSC in eight brain systems and mean age. This result can be reproduced using AAL1 (Extended Data Fig. 2-1) and Gordon's atlases (Extended Data Fig. 2-2). IVSC, individual variability in structural connectivity; SM, somatomotor network; VN, visual network; DMN, default mode network; DAN, dorsal attention network; Sub, subcortical system; FPN, frontoparietal network; VAN, ventral attention network; Lim, limbic network.
Figure 2-1
Quantify the IVSC map based on AAL atlas in 8 brain systems. A, Mean IVSC in 8 brain systems. B-I, The relationship between mean IVSC in 8 brain systems and mean age. AAL, automated anatomical labelling; IVSC, individual variability in structural connectivity; SM, somatomotor network; VN, visual network; DMN, default mode network; DAN, dorsal attention network; Sub, subcortical system; FPN, frontoparietal network; VAN, ventral attention network; Lim, limbic network. Download Figure 2-1, TIF file.
Figure 2-2
Quantify the IVSC map based on Gordon’s atlas in 7 brain systems. A, Mean IVSC in 8 brain systems. B-H, The relationship between mean IVSC in 8 brain systems and mean age. IVSC, individual variability in structural connectivity; SM, somatomotor network; VN, visual network; DMN, default mode network; DAN, dorsal attention network; FPN, frontoparietal network; VAN, ventral attention network; Lim, limbic network. Download Figure 2-2, TIF file.
To assess the robustness of the spatial distribution of IVSC, we performed a series of reproducible analyses. Whether to use spherical deconvolution-informed filtering of tractograms, and different weighting methods, had a minimal impact on the IVSC pattern. The IVSC map without spherical deconvolution-informed filtering of tractograms is shown in Extended Data Figure 1-4A, and the IVSC map based on networks weighted by the density of streamlines is shown in Extended Data Figure 1-4B. Their Spearman's r values with the IVSC map in the main results are 0.99 (p < 0.001) and 0.94 (p < 0.001), respectively.
Additionally, we assessed the influence of edge threshold selection. IVSC maps were calculated on the combined dataset after removing weak edges under different edge thresholds. The IVSC maps remained highly similar to the IVSC map without removing edges (Extended Data Fig. 1-5; all Spearman's r(246) values >0.96).
We also replicated the above analyses with different atlases: AAL1 atlas (Tzourio-Mazoyer et al., 2002) and a cortical parcellation created by Gordon and his colleagues (Gordon et al., 2016). The IVSC maps under different atlases were similar, with all maps showing higher variability in the medial prefrontal cortex, medial temporal cortex, anterior cingulate cortex, and middle frontal cortex and lower variability in the superior temporal cortex, the precuneus, the posterior cingulate cortex, and the regions near the central sulcus (Extended Data Fig. 1-6A,E). In addition, the relationship between IVSC and age was consistent across different atlases (Extended Data Fig. 1-6B–D,F–H). We mapped the AAL1 atlas into seven cortical networks and a subcortical system. Gordon's atlas was mapped into seven cortical networks only, because it does not include subcortical regions. The networks with the highest or lowest variability were consistent across different atlases, but the order of networks with moderate variability was slightly different across the three atlases (Extended Data Figs. 2-1A, 2-2A). The significant increases of the mean IVSC in all brain systems with age were reproduced using the AAL and Gordon's atlases (Extended Data Figs. 2-1B–I, 2-2B–H).
Association between IVSC and other properties of brain organization
Having identified the specific spatiotemporal pattern of IVSC, we next sought to understand the biological mechanism underlying this specific spatial distribution of IVSC. The human cortex, including the neocortex and allocortex, consists of several layers; each layer has a characteristic distribution of different neurons and connections with other cortical and subcortical regions. A previous study in mammals demonstrated that the existence of connections is not only related to the physical distance but is also related to the cytology of cortical areas (Goulas et al., 2019). In addition, myelin proteins can inhibit axon sprouting and stop the formation of structural connections by collapsing the tips of growing axons (Fields, 2008). Hence, we directly linked IVSC to these fundamental properties of brain organization. Using spatial permutation testing, we found that high IVSC was associated with low laminar differentiation (Fig. 3D; Spearman's r(210) = −0.46; p = 0.002), low myelin content (Fig. 3E; Spearman's r(210) = −0.44; p < 0.001), and low strength of short SC (Fig. 3F; Spearman's r(246) = −0.42; p < 0.001). To evaluate the stability of the correlation between IVSC and the strength of short SC, we filtered WM networks with different edge thresholds. Then, we calculated the IVSC and performed a correlation analysis for each edge threshold. We found that the correlation remained significant under different edge thresholds (Extended Data Fig. 3-1; −0.43 ≤ Spearman’s r(246) ≤ −0.48).
IVSC relates to fundamental properties of brain organization. A, Spatial distribution of the laminar differentiation identified by Paquola et al. (2019). B, Spatial distribution of myelin content measured by T1w/T2w mapping (Glasser and Essen, 2011). C, Spatial distribution of short connectivity strength. D–F, IVSC correlates with the laminar differentiation, myelin content, and strength of short connectivity. The correlation between IVSC and short connectivity strength was validated under different edge thresholds (Extended Data Fig. 3-1). IVSC, individual variability in structural connectivity; SC, structural connectivity.
Figure 3-1
The effect of edge thresholds on the correlation between individual variability in SC and strength of short SC across brain regions. To evaluate whether the correlation was stable across different edge thresholds. We recalculated individual variability in SC after filtering the weak edges whose weights were smaller than a threshold and then examined the correlation between individual variability in SC and the strength of short SC. We found that the relationship was stable across different edge thresholds. SC, structural connectivity. Download Figure 3-1, TIF file.
Association between IVSC and gene expression
Next, we explored the genetic mechanisms of IVSC using a public transcriptional dataset. For the AHBA, after controlling the quality of gene expression and filtering nonbrain-specific genes, the expression levels of 2,112 brain-specific genes across 123 left-hemispheric brain regions were obtained. Using PLS analysis, a multivariate statistical method to extract covarying components from matrices, we found a component of gene expression positively correlated with IVSC (Fig. 4B; Spearman's r(123) = 0.48; p = 0.005) and accounting for 26.1% of the spatial variance in IVSC. The component of gene expression was the sum of weighted gene expression, which we referred to as a gene score. The medial–temporal lobe, anterior cingulate cortex, and insula showed the highest gene scores, while the precentral gyrus, postcentral gyrus, and occipital gyrus showed the lowest gene scores (Fig. 4A).
IVSC relates to gene expression profiles. A, Spatial distribution of the gene expression score identified by PLS. B, IVSC was correlated with the gene expression score. C, The proportion of genes in a positive gene set preferentially expressed in seven different cell types. D, The proportion of genes in a negative gene set preferentially expressed in seven different cell types. E, Gene scores through human development. F–J, IVSC correlates with gene score in the fetus, infant, child, adolescent, and adult groups. IVSC, individual variability in structural connectivity; Astro, astrocyte; Micro, microglia; OPC, oligodendrocyte precursor; Oligo, oligodendrocyte; Endo, endothelial; Neuro-ex, excitatory neurons; Neuro-in, inhibitory neurons.
The relationship between IVSC and gene score was determined by the genes that contributed most to the gene score. The weights of genes reflected the shared variance between levels of gene expression and gene score. We defined the genes with the top 50% positive/negative weights as strongly contributing genes. To explore the biological processes of the positive and negative gene sets, we performed gene set enrichment analyses in both gene sets, respectively. We identified two biological processes, “positive regulation of AMPA glutamate receptor clustering” (p = 0.049) and “positive regulation of neuron projection arborization” (p = 0.049), in the positive gene set but revealed no biological process in the negative gene set.
Having found the biological process enriched in the positive gene set, we next investigated whether the genes associated with IVSC were preferentially expressed in specific cell types. Specifically, the proportion of genes that were preferentially expressed in each cell type was calculated. Permutation tests were used to assess statistical significance. Positive genes were significantly more expressed in astrocytes (p < 0.001) and oligodendrocyte precursors (p = 0.004; Fig. 4C), while negative genes were significantly more expressed in excitatory (p < 0.001) and inhibitory (p < 0.001) neurons (Fig. 4D).
In addition, we used the BrainSpan dataset to evaluate the change in the gene score throughout human development. The gene score of each tissue sample was calculated based on the gene weights derived from the PLS analysis. The tissue samples were separated into five life stages (fetus, infant, child, adolescent, and adult) in terms of age. The gene scores of tissue samples from the same life stage and same region were averaged to yield a gene score per region and life stage (Fig. 4E). The gene scores in most regions increased with development. We also used the BrainSpan dataset to evaluate the reproducibility of the association between IVSC and gene scores. We found that IVSC significantly correlated with gene scores in each of the life stages (Fig. 4F–J; fetus, Spearman's r(14) = 0.80; p = 0.003; infant, Spearman's r(14) = 0.63; p = 0.017; child, Spearman's r(14) = 0.55; p = 0.039; adolescent, Spearman's r(14) = 0.68; p = 0.021; and adult, Spearman's r(14) = 0.81; p = 0.002). Notably, correlation analyses were performed across 14 brain regions with gene expression levels.
Cognitive prediction with SC link between IVSC and individual difference in cognition
After identifying the spatiotemporal pattern of IVSC and uncovering its biological mechanism, we aimed to understand the impact of this specific IVSC pattern on cognition and investigate the relationship between IVSC and cognitive variation. To achieve this, we performed cognitive predictions using SC and examined the relationship between the selected features and IVSC. SC can be used to accurately predict executive function (Fig. 5A; Fold 1, MAE = 37.2; p < 0.001; Pearson's r(434) = 0.54; p < 0.001; Fold 2, MAE = 34.5; p < 0.001; Pearson's r(434) = 0.56; p < 0.001) and attention (Fig. 5D; Fold 1, MAE = 10.2; p < 0.001; Pearson's r(437) = 0.66; p < 0.001; Fold 2, MAE = 10.1; p < 0.001; Pearson's r(436) = 0.66; p < 0.001) in matched split-half samples with modified ridge regression. Repeated random 2F-CV produced similar results when predicting executive function (Fig. 5B; mean MAE = 37.9; mean Pearson's r(434) = 0.52) and attention (Fig. 5E; mean MAE = 9.8; mean Pearson's r(436) = 0.67) and thus suggested that the matched split of the data was unbiased. We also compared the performance of modified and traditional ridge regression models derived from repeated random 2F-CV. The modified ridge regression models outperformed the traditional ridge regression models in both tasks (Fig. 5B; executive function, two-sample t test, T(100) = −64.1; p < 0.001; Fig. 5E; attention function, two-sample t test, T(100) = −12.7; p < 0.001).
Prediction of executive function and attention benefits from individualized feature selection. To predict cognitive function, we used a modified ridge regression. Based on the model, we categorized SC into two kinds of edges: consistent edges and individualized edges (Extended Data Fig. 5-1). A, SC can be used to predict executive function with modified ridge regression. B, Modified RR outperformed RR in the task of predicting executive function. C, The percentage of consistent edges and individualized edges related to IVSC across brain regions when predicting executive function. D, SC can be used to predict attention function with modified ridge regression. E, Modified RR outperformed RR in predicting attention function. F, The percentage of consistent edges and individualized edges related to IVSC across brain regions when predicting attention function. MAE, mean absolute error; CV, cross-validation; RR, ridge regression; IVSC, individual variability in structural connectivity.
Figure 5-1
Model of modified ridge regression. A, Ridge regression was modified by integrating the ReLU activation function into it. B, For each edge, the number of participants who selected it as a feature was stated. Then, the edges that were selected by at least two-thirds of the participants were defined as consistent edges, and the edges that were selected by less than a third of participants were defined as individualized edges. Finally, we calculated the percentage of consistent edges and individualized edges for each node. ReLU, rectified linear unit. Download Figure 5-1, TIF file.
When predicting executive function and attention, the consistent SC is more than individualized SC, which means that the SC selected by models is majority common across individuals; only minority of SC is individual-specific (Fig. 5C,F). We found that the percentage of consistent edges was negatively correlated with IVSC, while the percentage of individualized edges was positively correlated with IVSC for predicting executive function (Fig. 5C; consistent edge, Pearson's r(246) = −0.14; p = 0.025; individualized edge, Pearson's r(246) = 0.18; p = 0.005) and attention (Fig. 5F; consistent edge, Pearson's r(246) = −0.23; p < 0.001; individualized edge, Pearson's r(246) = 0.31; p < 0.001).
Discussion
In this study, we used a human connectome method and multicenter data to delineate the spatiotemporal pattern of IVSC and explore its biological mechanisms. Building upon findings from this study, we identified a robust spatial distribution of IVSC: the limbic network shows the highest variability, association networks show moderate variability, and unimodal networks show the lowest variability. IVSC increases with increasing age in almost all cerebral areas. In addition, we provided evidence that IVSC is strongly linked to the microscopic anatomical properties of the brain and gene expression. Finally, we demonstrated that individualized feature selection improves the prediction. Taken together, these results offer a new perspective on brain organization and highlight the importance of individual effects in cognitive neuroscience research.
Spatiotemporal pattern of IVSC
The spatial distribution of IVSC we observed was moderately different from those in previous studies; however, consistent with the IVSC patterns presented in all these studies, we observed the lowest variability in the visual and somatomotor networks (Chamberland et al., 2017; Mansour et al., 2021; Sun et al., 2022). The inconsistent findings may mainly be due to selection bias. The sample size in one previous study was only nine, which is too small to represent the population. In addition, decreased IVFC was associated with increased genetic sharing magnitude (Gao et al., 2014). Therefore, including twins in the other two studies may lead to an underestimation of IVSC, especially in regions with high heritability. The interindividual variability in SC is lower than the variability in FC in adults; it exhibits a different spatial pattern mainly in the limbic network (i.e., temporopolar and orbitofrontal areas) that shows the highest variability in SC but moderately lower variability in FC. Intuitively, this disassociation is possible because resting-state functional connectivity not only reflects the underlying anatomical circuitry but also recapitulates the history of experience-driven coactivation (Lewis et al., 2009; Tambini et al., 2010).
Notably, the spatial distribution of IVSC is highly consistent with the sensory-fugal gradient, where one pole is occupied by the limbic areas and the opposite pole is the primary sensory–motor areas. The sensory–motor areas are polarized toward the extrapersonal environment and provide a gateway for the entry of sensory information into the cerebral neural circuitry and for outputting motor programs to bulbar and spinal motoneurons. In addition, the limbic areas are polarized toward the internal milieu and responsible for the regulation of memory, motivation, emotion, and autonomic–endocrine function (Mesulam, 1998). This distribution was proposed based on anatomical, physiological, and behavioral experiments in macaque monkeys and was demonstrated in the human brain with cytoarchitectonics, functional imaging, electrophysiological recordings, and the behavioral effects of focal lesions (Mesulam, 1998). The areas further away from the external input have more variable WM connections across human subjects. This characteristic of human networks may suggest that the flexibility of the “deep” cerebral layers, namely, trans-modal areas, supports the ability to adapt to changes in the environment, while the stability of sensory–motor areas plays a critical role in the accurate registration of new inputs.
Biological mechanisms of IVSC
Having identified the spatial distribution of IVSC, we sought to understand whether the macroscopic pattern of SC is constrained by the microscopic properties of brain structure. We observed that higher IVSC is associated with higher laminar differentiation and less myelin content. Recent evidence has demonstrated that laminar differentiation is related to the cortical plasticity–stability continuum of areas. The areas with lower laminar differentiation have more markers of synaptic plasticity (i.e., glial fibrillary acidic protein) and fewer markers of stability (i.e., parvalbumin), while the areas with the clearer laminar differentiation show the opposite trend (García-Cabezas et al., 2017). In addition, intracortical myelin is another marker that limits synaptic plasticity because the proteins in myelin, including Nogo-A (GrandPré et al., 2000; M. S. Chen et al., 2000), MAG (McKerracher et al., 1994), and OMgp (K. C. Wang et al., 2002; Huang et al., 2005), can lead to the collapse of tips of growing axons and thus prevent projections toward their target. Areas with higher plasticity are more easily influenced by life experience and may be more variable across populations. Hence, the heterogeneous spatial distribution of IVSC may be a macroscopic reflection of the differences in synaptic and axonal plasticity across different brain regions.
The relationship between IVSC and the strength of short connectivity suggests that the variation can be mainly sourced from the plastic change in weak short WM connectivity. However, although we used the most optimal processing pipeline and evaluated the correlation under different edge thresholds, this type of weak short connectivity may reflect noise inherent to the limitation of the tractography algorithm. Further studies are needed to prove this hypothesis. Importantly, we show that IVSC increases with aging, which is consistent with the findings in a study of functional connectivity (Ma et al., 2021). The age-related increase in IVSC indicates that the plastic change may last across the lifespan. Although we included only relatively young and healthy participants, the potential neuropathology of neurodegenerative disease may contribute to the part of the increase in variation.
Genetic mechanisms of IVSC
To further understand the genetic mechanism underlying the distribution of IVSC, we used expression data from AHBA and found two enriched biological processes: “positive regulation of AMPA glutamate receptor clustering” and “positive regulation of neuron projection arborization.” The AMPA receptor, a target of glutamate, has been proven to be important for modifying synaptic efficacy (Song and Huganir, 2002). Therefore, both biological processes are related to changes in neural pathways, which are the anatomical infrastructure for supporting learning. In addition, we observed that areas with higher IVSC are enriched for genetic signals of astrocytes and oligodendrocyte precursors. Interestingly, axonal electrical activity not only mediates the proliferation of oligodendrocyte precursors by controlling the production and release of growth factors (Barres and Raff, 1993) but also triggers the release of cytokine leukemia inhibitory factor by astrocytes, which promotes myelination (Ishibashi et al., 2006). Combining this evidence, the areas with higher IVSC may exhibit higher expression of genes involved in the formation of new pathways and the strengthening of existing pathways, indicating that these areas possess high plasticity. This conclusion is consistent with our findings from the correlation between IVSC and other microscopic properties. Although we used the data from different omics, the results suggest similar conclusions. In summary, heterogeneous spatial distribution of IVSC may be a macroscopic reflection of the regional differences in synaptic and axonal plasticity, likely resulting from the selective expression of the related genes.
Link between IVSC and variation in cognition
The modified model we proposed significantly improved the accuracy of predicting cognition. One possible reason for this improvement is that individuals recruit different neural circuits to perform cognitive tasks, especially advanced cognitive functions such as executive and attention functions. The modified model can select SC at the individual level to predict cognition, aligning with the hypothesis of individual variability in the neural circuitry of executive and attention function. Additionally, our findings suggest that variation in cognition may not only arise from differences in the quantity of brain measures within common neural circuits but may also be explained by the individual configuration of these neural circuits. fMRI studies have demonstrated a similar conclusion. For example, one study reported significant interindividual differences in brain functional parcellation, which is related to cognitive variation (Cui et al., 2020). Another study showed that higher intelligence is associated with the resting-state functional network that more closely resembles the task-related functional network (Schultz and Cole, 2016). Our results provide the complementary evidence from the perspective of SC for the link between cognition and neural configuration. The association between IVSC and percentage of individualized edges suggests that the contribution of individualized edges to cognitive prediction is not statistical noise but reflects underlying biological mechanisms. This strengthens our conclusion that individuals recruit different neural circuits to perform cognitive tasks. However, the improvement in cognitive prediction may simply be due to increased model complexity rather than the ability to select features at the individual level. Therefore, this conclusion needs to be verified in future studies.
The correlation between consistent/individualized edges and IVSC suggests that the distinct parts of neural circuits for a specific cognition are the edges connected to regions with high axonal and synaptic plasticity. Therefore, the specific IVSC pattern may partly result from the development of advanced cognition, during which individuals develop unique cognitive neural circuits, leading to the formation of different structural networks. Furthermore, the frontoparietal and ventral attention networks, which are involved in executive and attention functions, show both high IVSC and more individualized SC in cognitive prediction. This further supports that neural circuitry of executive and attention functions varies among the population. Collectively, our findings highlight the importance of brain mapping of cognition at the individual level for understanding neural mechanisms underlying cognitive function.
Potential practical applications of IVSC
Currently, precision medicine referring to tailoring treatment for individual is gradually becoming the future direction in healthcare. Understanding the spatial distribution of IVSC is likely to be crucial. For example, brain stimulation targets in areas with low IVSC may be suitable for group treatments, while targets in areas with high IVSC should be determined based on a patient's specific connectivity pattern. Furthermore, knowing the distribution of IVSC is likely beneficial for interpreting statistical maps. Traditional neuroimaging studies are more likely to obtain a significant result in areas with low variability, but not the areas with high variability. Therefore, the IVSC map might be used to correct this bias. Notably, these potential applications need to be systematically verified in future studies.
Limitation
Several methodological considerations in the present study should be noted. First, the structural networks we used to estimate IVSC were constructed based on a population-based atlas. However, the spatial topological properties of cortical regions vary across individuals (D. Wang et al., 2015; Kong et al., 2019; Han et al., 2020). Therefore, using a population-based atlas may be disadvantageous for estimating IVSC. Second, we regressed out intraindividual variability from interindividual variability, but the intraindividual variability and interindividual variability were estimated using datasets, which may also influence the estimation of IVSC. Third, the transcriptional data that we used were based on small samples of postmortem brains. More comprehensive microarray gene expression datasets are necessary for future studies.
Conclusion
In summary, we provide evidence that IVSC is constrained by microanatomical properties and gene expression. These findings may indicate that the pattern of IVSC is the reflection of regional differences in plasticity and the result of selective expression of related genes. Finally, our results also highlight the importance of individual connectivity patterns when explaining individual differences in cognition.
Footnotes
We thank all the volunteers for their participation in the study. This work was supported by the STI2030-Major Projects (2022ZD0213300, 2021ZD0200500), National Natural Science Foundation of China (32271145, 81871425, 81971585, 82271953, 82301608), Fundamental Research Funds for the Central Universities (2017XTCX04), the Open Research Fund of the State Key Laboratory of Cognitive Neuroscience and Learning (CNLZD2101, CNLYB2001), Key Research and Development Program of Hebei Provincial Department of Science and Technology (223777112D), Science and Technology Research and Development Plan of Chengde (202109A057), and Hebei Provincial Government-funded outstanding talent project.
The authors declare no competing financial interests.
- Correspondence should be addressed to Ni Shu at nshu{at}bnu.edu.cn.