Introduction

In neural circuitry, long-range (extracortical) interconnections among local (intracortical) microcircuits shape and constrain the large-scale functional organization of neural activity across the cortex[1–5]. The coupling of structural connectome (SC) and functional connectome (FC) varies greatly across different cortical regions reflecting anatomical and functional hierarchies[1, 6–9] and is regulated in part by genes[6, 8], as well as its individual differences relates to cognitive function[8, 9]. Despite its fundamental importance, our understanding of the changes in SC-FC coupling with development is currently limited. Specifically, the alterations in SC-FC coupling during development, its association with cognitive functions, and the underlying spatial transcriptomic mechanisms remain largely unknown.

Network modelling of the brain enables the characterization of complex information interactions at a system level and provides natural correspondences between structure and function in the cortex[7, 10]. Advances in diffusion MRI (dMRI) and tractography techniques have allowed the in vivo mapping of the white matter (WM) connectome (WMC), which depicts extracortical excitatory projections between regions[11]. The T1- to T2-weighted (T1w/T2w) ratio of MRI has been proposed as a means of quantifying microstructure profile covariance (MPC), which reflects a simplified recapitulation in cellular changes across intracortical laminar structure[6, 12–15]. Resting state functional MRI (rs-fMRI) can be used to derive the FC, which captures the synchronization of neural activity[16]. A variety of statistical[6, 8, 9], communication[1, 7], and biophysical[5, 17] models have been proposed to study the SC-FC coupling. The communication model is particularly useful because it not only depicts indirect information transmission but also takes into account biodynamic information within acceptable computational complexity[7, 18]. However, most studies have relied on WMC-derived extracortical communications as SC to predict FC, while ignoring the intracortical microcircuits, the MPC. In the present study, we propose that incorporating both intracortical and extracortical SC provides a more comprehensive perspective for characterizing the development of SC-FC coupling.

Previous studies in adults have revealed that the SC-FC coupling is strongest in sensory cortex regions and weakest in association cortex regions, following the general functional and cytoarchitectonic hierarchies of cortical organization[1]. This organization may occur due to structural constraints, wherein cortical areas with lower myelination and weaker WM connectivity tend to have more dynamic and complex functional connectivity[1, 8]. Large-scale association networks emerged over evolution by breaking away from the rigid developmental programming found in lower order sensory systems[19], facilitating regional and individual specialization[20]. In terms of developmental changes in SC-FC coupling, a statistical model-based study[9] identified positive age-related changes in some regions, while fewer regions exhibited negative changes. Furthermore, there is evidence that SC-FC coupling is linked to cognitive functions in healthy children[21], adults[8, 22] and patients[23], suggesting that it may be a critical brain indicator that encodes individual cognitive differences. Nonetheless, a more comprehensive investigation is needed to understand the precise pattern of SC-FC coupling over development and its association with cognitive functions.

Cortical SC-FC coupling is highly heritable[8] and related to heritable connectivity profiles[6], suggesting that the development of coupling may be genetically regulated. The Allen Human Brain Atlas (AHBA)[24] is a valuable resource for identifying genes that co-vary with brain imaging phenotypes and for exploring potential functional pathways and cellular processes via enrichment analyses[25–27]. For instance, a myeloarchitectural study showed that enhanced myelin thickness in mid-to-deeper layers is specifically associated with the gene expression of oligodendrocytes[28]. Another functional study found that the expression levels of genes involved in calcium ion-regulated exocytosis and synaptic transmission are associated with the development of a differentiation gradient[29]. However, the transcriptomic architecture underlying the development of SC-FC coupling remains largely unknown.

In this study, we analysed data obtained from the Lifespan Human Connectome Project Development (HCP-D)[30], which enrolled healthy participants ranging in age from 5.7 to 21.9 years. Our main objective was to investigate the SC-FC coupling of brain connectome and characterize its developmental landscapes. Specifically, we aimed to determine whether the SC-FC coupling encodes individual differences in cognition during development. Finally, we explored the genetic and cellular mechanisms underlying the development of SC-FC coupling of brain connectome. To assess the reproducibility of our findings, sensitivity and replication analyses were performed with a different tractography algorithm and a split-half independent validation method.

Results

We selected 439 participants (5.7 – 21.9 years of age, 207 males) in the HCP-D dataset who met our inclusion criteria: available high-quality T1/T2, dMRI, and rs-fMRI data that met the quality control thresholds. For each participant, we generated multiple connectomes using 210 cortical regions from the Human Brainnetome Atlas (BNA)[31], which comprised MPC, WMC, and FC. Intracortical connectivity was represented by MPC. According to the WMC, twenty-seven weighted communication models[7] were calculated to characterize geometric, topological, or dynamic connectivity properties. Further details on these models can be found in Text S1. After analysis, we found that communicability[32], mean first passage times of random walkers[33], and flow graphs (timescales=1) provided the optimal combination of extracortical connectivity properties because of significantly predicting FC (p<0.05, 1,000 spin test permutations, Table S1). We used these three models to represent the extracortical connectivity properties in subsequent discovery and reproducibility analyses (Figure S1).

Spatial pattern of cortical SC-FC coupling

We used SCs (MPC and three WMC communication models) to predict FC per node based on a multilinear model[1] (Figure 1), and quantified the nodewise SC-FC coupling as an adjusted coefficient of determination r² . We observed that the grouped SC-FC coupling varied across cortical regions (mean adjusted r² = 0.14 ± 0.08, adjusted r² range = [0.03, 0.45], Figure 2A), and regions with significant coupling were located in the middle frontal gyrus, precentral gyrus, paracentral lobule, superior temporal gyrus, superior parietal lobule, postcentral gyrus, cingulate gyrus, and occipital lobe (p<0.05, 1,000 spin test permutations, Figure 2B). Similar heterogeneous patterns of coupling were observed when categorizing cortical regions into seven functional subnetworks[34] (visual, somatomotor, dorsal attention, ventral attention, limbic, frontoparietal and default mode networks). In the visual, somatomotor, default mode and ventral attention networks, SC significantly predict FC variance (p<0.05, 1,000 spin test permutations, Figure 2C). The visual and somatomotor networks had higher coupling values than the other networks (p<0.05, Kruskal-Wallis ANOVA, Figure 2C). We further investigated the alignment between SC-FC coupling and three fundamental properties of brain organization: evolution expansion[35], myelin content[36], and functional principal gradient[37].

Figure 1.

Figure 2.

Our findings reveal a negative association between regional distribution of SC-FC coupling and evolution expansion (Spearman’s r=-0.52, p<0.001, 1,000 spin test permutations, Figure 2D), as well as with the functional principal gradient (Spearman’s r=-0.46, p<0.001, 1,000 spin test permutations, Figure 2F). Conversely, nodes exhibiting higher SC-FC coupling tended to exhibit higher myelin content (Spearman’s r=0.49, p<0.001, 1,000 spin test permutations, Figure 2E). In addition, the coupling pattern based on other models (using only MPC or only SCs to predict FC) and the comparison between the models were shown in Figure S2A-C.

Additionally, we applied Haufe’s inversion transform[38] to yield predictor weights of various SCs, where higher or lower values indicate stronger positive or negative correlations with FC. Our results demonstrated that different SCs had preferential contributions to FC variance across cortical regions to explain FC variance (p<0.05, FDR corrected, Kruskal-Wallis ANOVA, Figure 2G). Specifically, in the MPC, regions with positive correlation were the orbital gyrus, precentral gyrus, right middle temporal gyrus and temporoparietal junction, while regions with negative correlations were the left superior frontal gyrus, inferior parietal lobule and bilateral cingulate gyrus. Regarding WMC communication models, the communicability and flow graphs tended to stronger higher positive correlations in the visual, limbic and default mode networks, whereas the mean first passage time had stronger negative correlations in the somatomotor, limbic and frontoparietal networks.

Age-related changes in SC-FC coupling with development

To track changes in SC-FC coupling during development, we used a general linear model to assess the effect of age on nodal SC-FC coupling, while controlling for sex, intracranial volume, and in-scanner head motion. Our results revealed that the whole-cortex average coupling increased during development (β_age=1.05E-03, F=3.76, p=1.93E-04, r=0.20, p=3.20E-05, Figure 3A). Regionally, the SC-FC coupling of most cortical regions increased with age (p<0.05, FDR corrected, Figure 3B), particularly that in the frontal lobe, middle temporal gyrus, inferior temporal gyrus, parietal lobe, cingulate gyrus and lateral occipital cortex. Conversely, cortical regions with significantly decreased SC-FC coupling (p<0.05, FDR corrected, Figure 3B) were located in left orbital gyrus, left precentral gyrus, right superior and inferior temporal gyrus, left fusiform gyrus, left superior parietal lobule, left postcentral gyrus, insular gyrus, and cingulate gyrus. Age correlation coefficients distributed within functional subnetworks were shown in Figure 3C. Regarding mean SC-FC coupling within functional subnetworks, the somatomotor (β_age=2.39E-03, F=4.73, p=3.10E-06, r=0.25, p=1.67E-07, Figure 3E), dorsal attention (β_age =1.40E-03, F=4.63, p=4.86E-06, r=0.24, p=2.91E-07, Figure 3F), frontoparietal (β_age=2.11E-03, F=6.46, p=2.80E-10, r=0.33, p=1.64E-12, Figure 3I) and default mode (β_age=9.71E-04, F=2.90, p=3.94E-03, r=0.15, p=1.19E-03, Figure 3J) networks significantly increased with age and exhibited greater increase. No significant correlations were found between developmental changes in SC-FC coupling and the fundamental properties of cortical organization. Additionally, weights of different SCs varied with age, showing that MPC weight was positively correlated with age and that the weights of WMC communication models were stable (Figure S3-S6). The age-related patterns of SC-FC coupling based other coupling models were shown in Figure S2D-F.

Figure 3.

SC-FC coupling predicts individual differences in cognitive functions

As we found that SC-FC coupling can encode brain maturation, we next evaluated the implications of coupling for individual cognition using Elastic-Net algorithm[11]. After controlling for sex, intracranial volume and in-scanner head motion, we found the SC-FC coupling significantly predicted individual differences in fluid intelligence, crystal intelligence and general intelligence (Pearson’s r=0.3∼0.4, p<0.001, FDR corrected, Figure 4A). Furthermore, even after controlling for age, SC-FC coupling remained a significant predictor of general intelligence better than at chance (Pearson’s r=0.11± 0.04, p=0.01, FDR corrected, Figure 4A). For fluid intelligence and crystal intelligence, the predictive performances of SC-FC coupling were not better than at chance (Figure 4A). The predictive performances for other cognitive subscores are shown in Figure S7. To identify the regions with the greatest contributions to individual differences in age-adjusted general intelligence, we utilized Haufe’s inversion transform[38] to extract predictor weights across various regions. Our analysis revealed that SC-FC coupling within the prefrontal lobe, temporal lobe and lateral occipital lobe was the most predictive of individual differences in general intelligence (Figure 4B). In addition, we found that the weights of frontoparietal and default mode networks significantly contributed to the prediction of the general intelligence (p<0.01, 1,000 spin test permutations, Figure 4C).

Figure 4.

Transcriptomic and cellular architectures of SC-FC coupling development

We employed partial least square (PLS) analysis[39] to establish a link between the spatial pattern of SC-FC coupling development and gene transcriptomic profiles (Figure 5A) obtained from the AHBA using a recommended pipeline[40]. The gene expression score of the first PLS component (PLS1) explained the most spatial variance, at 22.26%. After correcting for spatial autocorrelation[41], we found a positive correlation (Pearson’s r=0.41, p=0.006, 10,000 spin test permutations, Figure 5B) between the PLS1 score of genes and the spatial pattern of SC-FC coupling development. In addition, we identified potential transcriptomic architectures using a Gene Ontology (GO) enrichment analysis of biological processes and pathway[42], analysing the significant positive and negative genes in PLS1. The positive weight genes (364 genes) were prominently enriched for “myelination”, “monoatomic cation transport”, “supramolecular fiber organization”, etc (p<0.05, FDR corrected, Figure 5C). The negative correlation genes (456 genes) were relatively weakly enriched in “cellular macromolecule biosynthetic process” and other pathways (p<0.05, FDR corrected, Figure 5C).

Figure 5.

To further investigate cell-specific expression patterns associated with SC-FC coupling development, the selected genes in the AHBA were agglomerated into seven canonical cell classes[43–48]: astrocytes, endothelial cells, excitatory neurons, inhibitory neurons, microglia, oligodendrocytes and oligodendrocyte precursors (OPCs). Our findings showed that the genes with positive weights were significantly expressed in oligodendrocytes (75 genes, p<0.001, permutation test, Figure 5D). The genes with negative weights were expressed in astrocytes (43 genes, p<0.001, permutation test, Figure 5D). Additionally, genes enriched in positive pathways were intensively overexpressed in oligodendrocytes, while genes enriched in three negative pathways were expressed in astrocytes, inhibitory neurons and microglia (p<0.05, permutation test, Figure S8).

Reproducibility analyses

Different parcellation templates

To evaluate the robustness of our findings to different parcellation templates, using the multimodal parcellation from the Human Connectome Project (HCPMMP)[49], we repeated the analyses of the cortical patterns of SC-FC coupling, correlation of age with SC-FC coupling, and gene weights. We observed a similar distribution in SC-FC coupling in which visual and somatomotor networks had higher coupling values than other networks (Figure S9A). The SC-FC coupling of most cortical regions increased with age (Figure S9B), and the significant regions were similar to those in the main findings (Figure S9C, p<0.05, FDR corrected). The gene weights of HCPMMP was consistent with that of BNA (r = 0.25, p < 0.001).

Different tractography strategies

To evaluate the sensitivity of our results to tractography strategies, we reconstructed fibres using deterministic tractography with a ball-and-stick model and generated a fibre number-weighted network for each participant. This same pipeline was employed for subsequent SC-FC coupling, prediction, and gene analyses. These two tractography strategies yielded similar findings, as indicated by significant correlations in the mean SC-FC coupling (r=0.85, p<0.001, spin test, Figure S10A), the correlation of between age and SC-FC coupling (r=0.79, p<0.001, spin test, Figure S10B), predictive weights on the general intelligence (r=0.85, p<0.001, spin test, Figure S10C), and gene weights (r=0.80, p<0.001, Figure S10D).

Split-half validation

To assess the reproducibility of our findings, we performed a split-half independent validation using the whole dataset (WD). Specifically, we randomly partitioned WD into two independent subsets (S1 and S2), and this process was repeated 1,000 times to mitigate any potential bias due to data partitioning. We then quantified SC-FC coupling, correlation between age and SC-FC coupling, and gene weights in S1 and S2 using the same procedures. Remarkably, we observed high levels of agreement among the datasets (S1, S2, and the WD) as demonstrated in Figure S11.

Discussion

In the present study, we characterized alterations of SC-FC coupling of brain connectome during development by combining intracortical and extracortical SC to predict FC based on the HCP-D dataset. We observed that SC-FC coupling was stronger in the visual and somatomotor networks than in other networks, and followed fundamental properties of cortical organization. With development, SC-FC coupling exhibited heterogeneous changes in cortical regions, with significant increases in the somatomotor, frontoparietal, dorsal attention and default mode networks. Furthermore, we found that SC-FC coupling can predict individual differences in general intelligence, mainly with the frontoparietal and default mode networks contributing higher weights. Finally, we demonstrated that the spatial heterogeneity of changes in SC-FC coupling with age was associated with transcriptomic architectures, with genes with positive weights enriched in oligodendrocyte-related pathways and genes with negative weights expressed in astrocytes.

Together, these findings characterized the spatial and temporal pattern of SC-FC coupling of brain connectome during development and the heterogeneity in the development of SC-FC coupling is associated with individual differences in intelligence and transcriptomic architecture.

Intracortical microcircuits are interconnected through extracortical WM connections, which give rise to richly patterned functional networks[1, 3]. Despite extensive research on this topic, the relationship between SC and FC remains unclear. Although many studies have attempted to directly correlate FC with the WMC, this correspondence is far from perfect due to the presence of polysynaptic (indirect) structural connections and circuit-level modulation of neural signals[2, 9, 16, 50]. Biological models can realistically generate these complex structural interconnections, but they have significant temporal and spatial complexity when solving for model parameters[51–54]. Communication models using the WMC integrate the advantages of different communication strategies and are easy to construct[18]. As there are numerous communication models, we identified an optimal combination consisting of three decentralized communication models based on predictive significance: communicability, mean first passage times of random walkers and flow graphs. We excluded a centralized model (shortest paths), which was not biologically plausible since it requires global knowledge of the shortest path structure[7, 55, 56]. In our study, we excluded the Euclidean distance and geodesic distance because spatial autocorrelation is inhibited. This study provides a complementary perspective (in addition to the role of WMC in shaping FC) that emphasizes the importance of intrinsic properties within intracortical circuit in shaping the large-scale functional organization of the human cortex. MPC can link intracortical circuits variance at specific cortical depths from a graph-theoretical perspective, enabling reflection of intracortical microcircuit differentiation at molecular, cellular, and laminar levels[6, 12–15]. Coupling models that incorporate these microarchitectural properties yield more accurate predictions of FC from SC[3, 57].

SC-FC coupling may reflect anatomical and functional hierarchies. SC-FC coupling in association areas, which have lower structural connectivity, was lower than that in sensory areas. This configuration effectively releases the association cortex from strong structural constraints imposed by early activity cascades, promoting higher cognitive functions that transcend simple sensori-motor exchanges[19]. A macroscale functional principal gradient[37, 58] in the human brain has been shown to align with anatomical hierarchies. Our study revealed a similar pattern, where SC-FC coupling was positively associated with evolutionary expansion and myelin content, and negatively associated with functional principal gradient during development. These findings are consistent with previous studies on WMC-FC coupling[9] and MPC-FC coupling[6]. Notably, we also found that the coupling pattern differed from that in adults, as illustrated by the moderate coupling of the sensorimotor network in the adult population[8]. SC-FC coupling is dynamic and changes throughout the lifespan[7], particularly during adolescence[6, 9], suggesting that perfect SC-FC coupling may require sufficient structural descriptors. Moreover, our results suggested that regional preferential contributions across different SCs lead to variations in the underlying communication process. Interestingly, the two extremes of regions in terms of MPC correlations corresponded to the two anchor points of the gradient[28]. The preferential regions in WM communication models were consistent with the adult results[7].

In addition, we observed developmental changes in SC-FC coupling dominated by a positive increase in cortical regions[9], broadly distributed across somatomotor, frontoparietal, dorsal attention, and default mode networks[9]. In a lifespan study, the global SC-FC coupling alterations with age were driven by reduced coupling in the sensorimotor network[7]. This finding is consistent across age ranges, indicating that sensorimotor coupling changes appear throughout development and ageing. Furthermore, we investigated the relationships of coupling alterations with evolutionary expansion and functional principal gradient but found no significant correlations, in contrast to a previous study[9]. These discrepancies likely arise from differences in coupling methods. We also found the SC-FC coupling with age across regions within subnetworks has more variability than the differences between networks, suggesting that the coupling with age is more likely region-dependent than network-dependent.

The neural circuits in the human brain support a wide repertoire of human behaviour[59]. Our study demonstrates that the degree of SC-FC coupling in cortical regions can significantly predict cognitive scores across various domains, suggesting that it serves as a sensitive indicator of brain maturity. Moreover, even after controlling for age effects, SC-FC coupling significantly predicted general intelligence, suggesting that it can partly explain individual differences in intelligence, as shown in previous studies[8]. In another study[9], positive correlations between executive function and SC-FC coupling were mainly observed in the rostro-lateral frontal and medial occipital regions, whereas negative associations were found in only the right primary motor cortex. While SC-FC coupling was not found to predict age-adjusted executive function in our study, we observed that the frontoparietal network and the default mode network specifically contributed higher positive prediction weights for general intelligence, whereas the somatomotor network had negative prediction weights[8]. The maturation of the frontoparietal network and default mode network continues into early adulthood, providing an extended window for the activity-dependent reconstruction of distributed neural circuits in the cross-modal association cortex[19]. As we observed increasing coupling in these networks, this may have contributed to the improvements in general intelligence, highlighting the flexible and integrated role of these networks.

Classic twin studies have reported that the heritability of coupling differs among cortical regions, with higher heritability in the visual network than in other cortical networks[8]. An inverse correlation between the pattern of SC-FC coupling and heritable connectivity profiles has been reported[6]. This led us to hypothesize that the development of SC-FC coupling may be influenced by the expression patterns of the genetic transcriptome across various cell types with different spatial distributions. Our findings suggest that the spatial development of SC-FC coupling is associated with underlying transcriptome structure. Specifically, genes positively associated with the development of SC-FC coupling were enriched in oligodendrocyte-related pathways. Oligodendrocytes, specialized glial cells in the central nervous system, play a crucial role in myelination by producing myelin sheaths that enable saltatory conduction and provide metabolic support to axons[60]. Defects in myelination have been linked to developmental disorders[61]. This seems to indicate that significant alterations in SC-FC coupling during development may reflect neural plasticity, such as activity-dependent myelination of axons connecting functionally coupled regions[62, 63]. Conversely, we found that genes negatively correlated with SC-FC coupling were enriched in two specific gene pathways within astrocytes, inhibitory neurons and microglia. Both astrocytes and microglia have been implicated in synaptic pruning, a critical developmental process for the formation of fully functional neuronal circuits that eliminates weak and inappropriate synapses[64–66]. Importantly, the precise establishment of synapses is crucial for establishing the intercellular connectivity patterns of GABAergic neurons[67]. These findings suggest that the subtle alterations observed in SC-FC coupling are closely associated with the refinement of mature neural circuits.

Several methodological issues must be addressed. First, we implemented a conservative quality control procedure to address head motion, which unavoidably resulted in the loss of some valuable data. Given the confounding influence of head motion in fMRI studies, especially those involving developing populations, we applied censoring of high-motion frames and included motion as a covariate in the GLM analysis and cognitive prediction to minimize its effects[7, 59, 68, 69]. Second, Second, although we observed SC-FC coupling across development by integrating intra- and extracortical SC to predict FC, it is worth noting that combining deep learning models[2], biophysical models[5, 17], or dynamic coupling[3, 13] perspectives may provide complementary insights. Third, the appropriateness of structurally defined regions for the functional analysis is also a topic of important debate. Fourth, we focused solely on cortico-cortical pathways, excluding subcortical nuclei from analysis. This decision stemmed from the difficulty of reconstructing the surface of subcortical regions[70] and characterizing their connections using MPC technique, as well as the challenge of accurately resolving the connections of small structures within subcortical regions using whole-brain diffusion imaging and tractography techniques[71, 72]. In addition, the reconstruction of short connections between hemispheres is a notable challenge. Fifth, it is important to acknowledge that changes in gene expression levels during development may introduce bias in the results. Finally, validation of sensitivity across independent datasets is a crucial step in ensuring the reliability of our results. To address this, we employed an alternative split-half validation strategy and the results supported the reliability of the current findings. However, future verification of current findings on independent datasets are still needed.

Conclusions

Overall, this study sheds light on the development of SC-FC coupling in the brain and its relationship to cognitive function and gene expression patterns. The results improve our understanding of the fundamental principles of brain development and provide a basis for future research in this area. Further investigations are needed to fully explore the clinical implications of SC-FC coupling for a range of developmental disorders.

Materials and Methods

Participants

We selected 439 participants (207 males, mean age = 14.8 ± 4.2 years, age range = [5.7, 21.9]) from the HCP-D Release 2.0 data (https://www.humanconnectome.org/study/hcp-lifespan-development) after conducting rigorous checks for data completeness and quality control.

The HCP-D dataset comprised 652 healthy participants who underwent multimodal MRI scans and cognitive assessments, and the detailed inclusion and exclusion criteria for this cohort have been described in[30]. All participants or their parents (for participants under the age of 18 years) provided written informed consent and assent. The study was approved by the Institutional Review Board of Washington University in St. Louis.

Cognitive scores

We included 11 cognitive scores which were assessed with the National Institutes of Health (NIH) Toolbox Cognition Battery (https://www.healthmeasures.net/exploremeasurement-systems/nih-toolbox), including episodic memory, executive function/cognitive flexibility, executive function/inhibition, language/reading decoding, processing speed, language/vocabulary comprehension, working memory, fluid intelligence composite score, crystal intelligence composite score, early child intelligence composite score and total intelligence composite score. Distributions of these cognitive scores and their relationship with age are illustrated in Figure S12.

Imaging acquisition

The MRI data were obtained with a Siemens 3T Prisma with a 32-channel phased array head coil, and detailed imaging parameters are available in[73]. High-resolution T1w images were acquired using a 3D multiecho MPRAGE sequence (0.8 mm isotropic voxels, repetition time (TR)/inversion time (TI) = 2500/1000 ms, echo time (TE) = 1.8/3.6/5.4/7.2 ms, flip angle = 8°, up to 30 reacquired TRs). The structural T2w images were collected with a variable-flip-angle turbo-spin‒echo 3D SPACE sequence (0.8 mm isotropic voxels, TR/TE = 3200/564 ms, up to 25 reacquired TRs). The dMRI scans included four consecutive runs with a 2D 4× multiband spin‒ echo echo-planar imaging (EPI) sequence (1.5 mm isotropic voxels, 185 diffusion directions with b = 1500/3000 s/mm² and 28 b = 0 s/mm² volumes, TR = 3.23 s, flip angle = 78°). The rs-fMR images were acquired using a 2D 8× multiband gradient-recalled echo EPI sequence (2.0 mm isotropic voxels, TR/TE = 800/37 ms, flip angle = 52°). Each rs-fMRI scan duration was 26 minutes (four runs of 6.5 minutes) for participants over 8 years old and 21 minutes (six runs of 3.5 minutes) for participants who were 5∼7 years old.

Imaging preprocessing

All structural, diffusion and functional images underwent minimal preprocessing[70]. We specifically processed dMRI data referring to the publicly available code from https://github.com/Washington-University/HCPpipelines since the HCP-D has not released preprocessed dMRI results. Briefly, structural T1w and T2w images went through gradient distortion correction, alignment, bias field correction, registration to Montreal Neurological Institute (MNI) space, white matter and pial surface reconstruction, segment structures, and surface registration and downsampling to 32k_fs_LR mesh. A T1w/T2w ratio image, which indicates intracortical myelin, was produced for each participant[36]. The BNA[31] was projected on native space according to the official scripts (http://www.brainnetome.org/resource/) and the native BNA was checked by visual inspection. Regarding fMRI data, the preprocessing pipeline included spatial distortion correction, motion correction, EPI distortion correction, registration to MNI space, intensity normalization, mapping volume time series to 32k_fs_LR mesh, and smoothing using a 2 mm average surface vertex. Following our previous methodological evaluation study[11], the dMRI procedures consisted of intensity normalization of the mean b0 image, correction of EPI distortion and eddy current, motion correction, gradient nonlinearity correction, and linear registration to T1w space.

Network computation

Microstructure profile covariance (MPC)

The MPC can capture cytoarchitectural similarity between cortical areas[12]. We first reconstructed 14 cortical surfaces from the white matter to the pial surface using a robust equivolumetric model[12, 74]. Then, the T1w/T2w ratio image was used to sample intracortical myelin intensities at these surfaces. We averaged the intensity profiles of vertices over 210 cortical regions according to the BNA[31]. Finally, we computed pairwise partial correlations between regional intensity profiles, while controlling for the average intensity profile. After removing negative correlations, we used Fisher’s r-to-z-transformation to generate an individual MPC.

White matter connectome (WMC)

Following our previous methodological evaluation study[11], the ball-and-stick model estimated from the bedpostx command-line in the FDT toolbox of FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT) was used to estimate fibre orientations (three fibres modelled per voxel)[75–78]. The BNA atlas was applied to individual volume space by inverse transformation derived from preprocessed steps. Next, probabilistic tractography (probtrackx)[77, 79] was implemented in the FDT toolbox to estimate the probability of connectivity between two regions by sampling 5,000 fibres for each voxel within each region, correcting for distance, dividing by the total fibres number in source region, and calculating the average bidirectional probability[11]. Notably, the connections in subcortical areas were removed. A consistency-based thresholding approach (weight of the coefficient of variation at the 75th percentile) was used to remove spurious connections, and retain consistently reconstructed connections across subjects[9, 80]. Finally, twenty-seven communication models[7] were subsequently derived from the weighted probabilistic network, including shortest path length (gamma values = {0.12, 0.25, 0.5, 1, 2, 4}), communicability[32], cosine similarity, search information (weight-to-cost transformations = {0.12, 0.25, 0.5, 1, 2, 4})[81], path transitivity (weight-to-cost transformations = {0.12, 0.25, 0.5, 1, 2, 4})[55], matching index[82], greedy navigation[83], mean first passage times of random walkers[33], and flow graphs (timescales = {1, 2.5, 5, 10})[84]; for more details see Text S1.

Functional network (FC)

To further clean the functional signal, we performed frame censoring, regressed out nuisance variables (including white matter, cerebrospinal fluid, global signal, and 12 motion parameters), and executed temporal bandpass filtering (0.01∼0.1 Hz). Specifically, we identified censored frames with motion greater than 0.15 mm[7] based on the Movement_RelativeRMS.txt file. We flagged one frame before and two frames after each censored frame, along with any uncensored segments of fewer than five contiguous frames, as censored frames as well[69]. We discarded fMRI runs with more than half of the frames flagged as censored frames, and excluded participants with fewer than 300 frames (less than 4 minutes). The nuisance variables were removed from time series based on general linear model. We averaged the time series of vertices into 210 cortical regions according to the BNA[31]. We then computed pairwise Pearson’s correlations between regional time series, and applied Fisher’s r-to-z-transformation to the resulting correlations to generate individual FC.

Quality control

The exclusion of participants in the whole multimodal data processing pipeline was depicted in Figure S13. In the context of fMRI data, we computed Pearson’s correlation between motion and age, as well as between the number of remaining frames and age, for the included participants aged 5 to 22 years and 8 to 22 years, respectively. These correlations were presented in Figure S14.

SC-FC coupling

A multilinear model[1] was constructed to examine the relationship of individual nodewise SC profiles and FC profiles. For a given node, the predictive variable was nodal SC S = {s₁, s₂, ⋯, s_i, ⋯, s_n}, s_i ∈ R^m where s_i is the ith SC profiles, n is the number of SC profiles, and m is the node number. The nodal functional profile f is the dependent variable.

where the intercept b₀ and regression coefficients b_i are estimated model parameters. For each participant, goodness of fit per node represents the nodal coupling between SC and FC, quantified as the adjusted coefficient of determination[7]

where R² is the unadjusted coefficient of determination, N_c is the number of connection (N_c = 245 for BNA), and N_p is the number of predictors.

In the present study, WMC communication models that represented diverse geometric, topological, or dynamic factors, were used to explain nodal FC variation. Notably, too many predictors will result in overfitting and blindly increase the explained variance. And covariance structure among the predictors may lead to unreliable predictor weights. Thus, we applied Haufe’s inversion transform[38] to address these issues and identified reliable communication mechanisms. Specifically, we used all twenty-seven communication models to predict FC at the node level for each participant. We applied Haufe’s inversion transform[38] to obtain predictor weights for each model, with higher or lower values indicating stronger positive or negative correlations with FC. Next, we generated 1,000 FC permutations through a spin test[85] for each nodal prediction in each subject and obtained random distributions of model weights. These weights were averaged over the group and were investigated the enrichment of the highest weights per region to assess whether the number of highest weights across communication models was significantly larger than that in a random discovery.

The significant communication models were used to represent WMC communication properties and to predict functional profiles in conjunction with MPC as structural profiles (predictors). To test the significance of the resulting adjusted R² values and system-specific of coupling, we generated a null predictive model using a spin test[85] with 1,000 spatially-constrained repetitions. We also used Kruskal-Wallis nonparametric one-way analysis of variance (Kruskal-Wallis ANOVA) to compare coupling differences between systems. To investigate the contributions of various structural predictors, we applied Kruskal-Wallis ANOVA to test the predictive weights derived by Haufe’s inversion transform, identifying optimal predictors across regions. We corrected for multiple comparisons using FDR correction. Additionally, we used a general linear model to explore age-related developmental patterns of SC-FC coupling, while controlling for sex, intracranial volume, and in-scanner head motion. Similarly, the system-specific significance of coupling alteration was calculated based on the 1,000 repetitions of the spin test. In addition, we have constructed the models using only MPC or SCs to predict FC, respectively. Spearman’s correlation was used to assess the consistency between spatial patterns based on different models.

We examined the associations of SC-FC coupling and its developmental pattern with evolution expansion[35], myelin content[36], and functional principal gradient[37]. Spearman’s correlation analyses were used to quantify the strength of correlations, with significance corrected for spatial autocorrelation with 1,000 repetitions of the spin test.

Prediction of cognitive function

Based on our predictive evaluation work[11], the Elastic-Net algorithm was applied to predict cognitive performance using nodal SC-FC coupling, which tends to yield robust prediction performance across various dimensions of cognitive tasks. The objective function is as follows:

where X = {x₁, x₂, …, x_n} represents an observation set (e.g., SC-FC coupling) with a sample size of n, and Y = {y₁, y₂, …, y_n} is a label set (e.g., cognitive measure). The model solves the fitting coefficient W = (w₁, w₂, …, w_m) under the minimization objective function L(Y, f(X, W)). The L1 regularization term |·| and L2 regularization term ||·||² constraint the fitting coefficient to ensure model generalization ability. α represents regularization strength, controlling the compression loss scale, and β denotes a trade-off parameter between the L1 and L2 terms.

We employed a nested fivefold cross validation (CV) framework comprising an external CV and an internal CV[11]. In the external CV, observations were randomly partitioned into five folds, with four of them included in the training set used to develop the model and the remaining fold used as a testing set to assess the predictive accuracy of the model. This process was repeated 100 times, and the final model performance was evaluated by averaging the predictive accuracy across the 100 models. In the internal CV, the hyperparameter spaces were first defined as α ∈ {x|x = 2ⁿ, n ∈ Z, n ∈ [−10,5]} and β ∈ {x|x = 0.1n, n ∈ Z, n ∈ [0,10]} . Then, the training set was further divided into five folds. Four folds composed the internal training set, which was used to generate models by successively applying 16×11 hyperparametric combinations, and the remaining fold was defined as the validation set and used to find the optimal combination. Subsequently, we retrained the model on the training set using the optimal hyperparametric combination and assessed its predictive performance on the testing set by performing Pearson’s correlation analyses of the relationship between the predicted and labelled values.

Prior to applying the nested fivefold cross-validation framework to each behaviour measure, we regressed out covariates including sex, intracranial volume, and in-scanner head motion from the behaviour measure[59, 69]. Specifically, we estimated the regression coefficients of the covariates using the training set and applied them to the testing set. This regression procedure was repeated for each fold. Additionally, we conducted control analyses using age-adjusted behavioral measures to investigate the effect of age on the predictive performance of SC-FC coupling.

To evaluate whether our model performed better than at chance on each behaviour measure, we performed 1,000 permutation tests by randomly shuffling the behaviour measure across participants, generating a null model of predicted performance using the same procedures. We then used the corrected resampled t test to determine statistical significance[86, 87]. We corrected for multiple comparisons using FDR correction. For model interpretability, we applied Haufe’s inversion transform[38] to obtain predicted weights for various brain regions. The significance of the weights for each system was assessed by comparing them to those generated by a spin test[85] with 1,000 repetitions.

Association between alterations of SC-FC coupling and gene expression

We preprocessed the anatomic and genomic information of the Allen Human Brain Atlas (AHBA) dataset following a recommended pipeline[40]. Specifically, we used FreeSurfer (https://surfer.nmr.mgh.harvard.edu/fswiki/) to generate preprocessed structural data for each donor and projected the BNA template onto native fsaverage space using official scripts (http://www.brainnetome.org/resource/). Finally, we produced an averaged gene expression profile for 10,027 genes covering 105 left cortical regions. Restricting analyses to the left hemisphere will minimize variability across regions (and hemispheres) in terms of the number of samples available[40].

PLS analysis[39] was performed to mine the linear association between the spatial development pattern of SC-FC coupling and gene expression profiles. We used absolute values of the correlation between age and SC-FC coupling in 105 left cortical regions as predicted variables and the gene expression profiles of the corresponding regions as predictor variables. Pearson’s correlation coefficient was calculated to determine the association between the PLS score and the absolute correlation value between age and SC-FC coupling. To correct for spatial autocorrelation, we compared the empirically observed value to spatially constrained null models generated by 10,000 spin permutations[85]. We then transformed the gene weight on PLS1 into a z score by dividing the standard deviation of the corresponding weights estimated from bootstrapping, and ranked all genes accordingly. We identified significant genes at a threshold of p < 0.05 and classified them as having positive or negative gene weights. To understand the functional significance of these genes, we performed gene functional enrichment analysis (GO analysis of biological processes and pathways) using Metascape[42]. We focused on the selected genes with positive or negative weights and retained enrichment pathways with an FDR-corrected < 0.05.

To investigate the cell type-specific expression of the selected genes, we assigned them to 58 cell types derived from five studies[43–47] focusing on single-cell research using the human postnatal cortex. To avoid potential bias in cell-type assignment, we grouped these cell types into seven canonical classes: astrocytes, endothelial cells, excitatory neurons, inhibitory neurons, microglia, oligodendrocytes, and OPCs[48, 88]. We generated a null model by performing 10,000 random resamplings of genes within each cell type. We then tested the significance of our results against this null model. Additionally, we subjected the genes associated with each enriched term to the same analysis to explore the specificity of the cell type.

Reproducibility analyses

To evaluate the robustness of our findings under different parcellation templates, we computed MPC, SCs (WMC, communicability[32], mean first passage times of random walkers[33], and flow graphs (timescales=1)) and FC using the multimodal parcellation from the Human Connectome Project (HCPMMP)[49]. We used the multilinear model to examine the association of individual nodewise SC profiles and FC profiles. Then, a general linear model was used to explore age-related developmental patterns of SC-FC coupling, while controlling for sex, intracranial volume, and in-scanner head motion. We corrected for multiple comparisons using FDR correlation. Finally, we produced an averaged gene expression profile for 10,027 genes covering 176 left cortical regions based on HCPMMP and obtained the gene weights by PLS analysis. We performed Pearson’s correlation analyses to assess the consistency of gene weights between HCPMMP and BNA.

To evaluate the sensitivity of our results to deterministic tractography, we used the Camino toolbox (http://camino.cs.ucl.ac.uk/) to reconstruct fibres with a ball-and-stick model estimated from bedpostx results[78] and to generate a fibre number-weighted network using the BNA atlas. We then calculated the communication properties of the WMC including communicability, mean first passage times of random walkers, and flow graphs (timescales=1). The same pipeline was used for subsequent SC-FC coupling, prediction, and gene analysis. To assess the consistency of our results between deterministic tractography and probabilistic tractography, we performed Pearson’s correlation analyses with significance corrected for spatial autocorrelation through 1,000 repetitions of the spin test.

To evaluate the generalizability of our findings, we adopted a split-half cross-validation strategy by randomly partitioning the whole dataset (WD) into two independent subsets (S1 and S2). This process was repeated 1,000 times to minimize bias due to data partitioning. Based on MPC, three communication properties of the WMC, and FC, we then used the same procedures to quantify SC-FC coupling, the correlation between age and SC-FC coupling and gene weights in both S1 and S2. Finally, we assessed the consistency of results by calculating Pearson’s correlation coefficients of the relationships between S1 and WD, S2 and WD, and S1 and S2.

Data and Code availability

The HCP-D 2.0 release data that support the findings of this study are publicly available at https://www.humanconnectome.org/study/hcp-lifespan-development.R4.1.2 software (https://www.r-project.org/) was used to construct the general linear model. MATLAB scripts used for preprocessing of the AHBA dataset can be found at https://github.com/BMHLab/AHBAprocessing. Python scripts used to perform PLS regression can be found at https://scikit-learn.org/stable/modules/generated/sklearn.cross_decomposition.PLSRegression.html#sklearn.cross_decomposition.PLSRegression. The minimal preprocessing pipelines can be accessed at https://github.com/Washington-University/HCPpipelines. The code relevant to this study can be accessed through the following GitHub repository: https://github.com/FelixFengCN/SC-FC-coupling-development.

Credit Authorship Contribution Statement

G.F. performed acquisition and analysis of data, contributed new analytic tools, drafted and revised the paper. Y.W., W.H., H.C. and J.C. performed acquisition and analysis of data. N.S. contributed to the design of the work, performed analysis of data, and revised the paper.

Acknowledgements

The authors thank all the volunteers for their participation in the study and anonymous reviewers for their insightful comments and suggestions. This work was supported by the STI2030-Major Projects (2021ZD0200500, 2022ZD0213300), National Natural Science Foundation of China (32271145, 81871425), Fundamental Research Funds for the Central Universities (2017XTCX04), Open Research Fund of the State Key Laboratory of Cognitive Neuroscience and Learning (CNLZD2101). Data in this publication were provide (in part) by the Human Connectome Project-Development (HCP-D), which is supported by the National Institute Of Mental Health of the National Institutes of Health under Award Number U01MH109589 and by funds provided by the McDonnell Center for Systems Neuroscience at Washington University in St. Louis.

Financial Disclosures

There are no conflicts of interest including any financial, personal, or other relationships with people or organizations for any of the authors related to the work described in the article.