Dual-axis myelination covariance drives the functional connectivity emergence during infancy

Abstract

The mechanisms linking structural maturation to the emergence of functional networks in the perinatal brain remain unresolved. While prevailing models attribute functional connectivity to white matter myelination, neonates paradoxically exhibit adult-like resting-state networks despite profoundly immature white matter tracts. Here, we proposed gray matter myelination covariance as a critical basis of early functional connectivity emergence. We introduced a dual-axis myelination covariance framework and derived a myelination-function coupling (MFC) index specific to the newborn brain. Results revealed that the MFC exhibited distinct spatial patterns dominated by primary sensory and motor cortices, increased with age, and showed a distance-dependent strength. Crucially, neonatal MFC patterns showed a strong spatial correlation with gene expression profiles implicated in neurovascular coupling and specifically predicted later behaviors. These findings suggest that during infancy, the integration of brain function is not initially dominated by only the white matter connections but is also shaped by the synchrony of intracortical microstructure that reflects shared developmental trajectories, which offers a framework for understanding the formation of the developmental connectome.

Fig. 1. The research pipeline of this study.

a Basic question of this study. b The brain imaging data and projection to the cortical surface. c The proposed workflow of myelination-function coupling (MFC) based on dual-axis covariance. d Developmental patterns of the MFC index. e Spatial distance dependence analysis of the MFC map. f Genetic decoding using transcriptional enrichment analysis. g The association analysis between behavior outcomes at 18-month-age and brain MFC at birth. MATLAB R2022b (MathWorks, Natick, MA, USA) and Connectome Workbench were used for visualizations.

Fig. 2. Spatial pattern of myelination-function coupling index.

a Linear regression between mean functional connectivity (mFC) and group-level myelination covariance (gMC) or subject-level myelination covariance (sMC) for the connections of intra-hemisphere or whole connectome (n = 18,438,448 for hemisphere; n = 36,881,166 for whole brain), with the background displaying a density distribution map. The exact p equals zero for each of the four associations. b The spatial pattern of the group-averaged MFC with dual-axis myelination covariance. c The spatial patterns of gMFC and sMFC. gMFC: the group-level MFC derived from dual-axis framework only using gMC; sMFC: the subject-level MFC derived from dual-axis framework only using sMC. d System-level distributions of the three coupling indices (displayed in colored histograms for each system), most of which are significantly different from the null distribution (displayed as grayscale histograms), and vertical lines in the graph represent the mean values. * p_spin < 0.05; ** p_spin < 0.01; *** p_spin < 0.001. The investigated systems are visual (VIS), sensorimotor (SM), frontoparietal (FP), default mode (DMN), limbic (LIM), dorsal attention (DA), and ventral attention (VA) systems. Source data are provided as a Source Data file. Analyses and visualization were performed using custom scripts in MATLAB R2022b (MathWorks, Natick, MA, USA) together with Connectome Workbench.

Fig. 3. Typical development of the myelination-function coupling.

a The spatial-temporal maps of MFC for a series of PMA windows. b Growth-rate coefficient of MFC during infancy derived from GAM analysis with p_FDR < 0.05. c Growth-rate coefficient of gMFC and sMFC with p_FDR < 0.05. d The distributions and comparisons of the growth-rate coefficient at significant vertices in GAM of the three coupling indices, and the statistical differences between them are all significant (* p_FDR < 0.05; ** p_FDR < 0.01; *** p_FDR < 0.001, the number of significant vertices for each system and each index is listed in Supplementary Table 4). Boxes display the interquartile range (IQR: lower hinge: 25th percentile; upper hinge: 75th percentile; center line: median). e Radar chart demonstrates the proportion of significantly developed areas for each system, with solid line for MFC, large dashed lines for gMFC and small dashed lines for sMFC. The investigated systems are visual (VIS), sensorimotor (SM), frontoparietal (FP), default mode (DMN), limbic (LIM), dorsal attention (DA), ventral attention (VA) systems. Source data are provided as a Source Data file. Analyses and visualization were performed using custom scripts in MATLAB R2022b (MathWorks, Natick, MA, USA) together with Connectome Workbench.

Fig. 4. Distance dependence of the MFC patterns and the growth rate.

a Spatial patterns of MFC from local connections to remote connections. b MFC curves (solid lines, with locally weighted scatterplot smoothing) and growth rate curves (dashed lines) against the distance for each system, with the dark colored curves and light colored curves respectively representing the right and left hemispheres respectively. The shadows for solid lines show the upper 60th percentile and the lower 40th percentile from n = 364. c Linear fitting coefficients at each vertex between the growth rate of MFC and cortical distance, positive values represent a higher growth rate for long-distance coupling, negative values represent a lower growth rate for short-distance coupling. d Averaged map of individual MFC with distance regressed. e Averaged map of individual MFC with local connections removed. Source data are provided as a Source Data file. Analyses and visualization were performed using custom scripts in MATLAB R2022b (MathWorks, Natick, MA, USA) together with Connectome Workbench.

Fig. 5. Birth effect on the MFC patterns.

a Regression performance of model II (GA and PNA combined) is greater than that of model I (only PMA) at the vertex level. b Distribution of GA coefficient and PNA coefficient at vertex level. c The regression performances ($R_{adj}^2$ of model I and model II in the left panels) or the coefficients (GA and PNA in the right panels) at significant fitting vertices with p_FWE < 0.05 for each system are illustrated in mean $\pm$SD. All group comparisons are significant (p_FDR < 0.05, two-sided) and t-values are illustrated by the size of gray circular bubble. The sample size which equals to the number of significant vertices for each bar are listed in Supplementary Table 7). d Age information of the preterm and full-term neonate subsets, including PMA, GA, and PNA, are showed as the mean $\pm$SE. GA and PNA show difference between the preterm (n = 83) and the age-matched full-term (n = 83) with (GA: t = -18.58, p = 3.4 $\times$ 10^-42, two sided; PNA: t = 14.46, p = 4.5 $\times$ 10⁻³¹, two sided). e Two-sample T-test indicates that the averaged vertex-level values of individual MFC of full-term group are higher than that of preterm group (t = 23.42, p = 1.8 $\times$ 10^-119, two-sided). Boxes display the interquartile range (IQR: lower hinge: 25th percentile; upper hinge: 75th percentile; center line: median). f Spatial maps of averaged vertex-level MFC corresponding to panel e. GA gestational age, PNA postnatal age, PMA postmenstrual age. Source data are provided as a Source Data file. Analyses and visualizations here were performed using MATLAB R2022a (MathWorks Inc., Natick, MA, USA).

Fig. 6. Transcriptomic decoding of MFC.

a Variance explained for the first 10 components derived from the PLS regression analysis. Permutation test is used by randomly shuffling the orders of regional MFC 5,000 times to assess the statistical significance of PLS1. b Pearson’s correlation between MFC score and gene expression score (n = 11, r = 0.91, p = 0.0001, CI = [0.90, 0.92], two-sided). Red line and shaded represent the linear regression line with 95% confidence intervals. c Ranked PLS1 loadings of significantly associated genes. d Enrichment analysis of PLS+ genes (orange bars, Top 20 GO enriched clusters) and PLS- genes (blue bars) associated with MFC, all -log(P-value) > 2.22. e The intra- and inter-cluster similarities of enriched terms of PLS+ genes are visualized by the Metascape enrichment network. CC Cellular Component, BP Biological Process, MF Molecular Function. Source data are provided as a Source Data file. Analyses and visualization were performed using custom scripts in MATLAB R2022b (MathWorks, Natick, MA, USA).

Fig. 7. Relationship between MFC and the behavior scores at 18 months of age.

a Significant association between behavioral and MFC composite scores under the PLS1 (n = 278, adjusted R² = 0.026, p_perm = 0.006, two-sided), with the behavioral and MFC loadings. Red line and shade represent the linear regression line with 95% confidence intervals. Permutation test is used by randomly shuffling the orders of regional MFC 5000 times to assess the association. b Significant association between behavioral and FCS composite scores in the first latent component (n = 278, adjusted R² = 0.047, p_perm = 0.001, two-sided), with the behavioral and FCS loadings. Red line and shade represent the linear regression line with 95% confidence intervals. Permutation test is used by randomly shuffling the orders of regional MFC 5,000 times to assess the association. c Significant association between the motor domain and MFC of motor network composite score in the first latent component (n = 278, adjusted R² = 0.045, p_perm = 0.041, two-sided), with the domain loading and the MFC loadings. Red line and shade represent the linear regression line with 95% confidence intervals. Permutation test is used by randomly shuffling the orders of regional MFC 5000 times to assess the association. d The monotonic relationship between the R² and the distance (n = 9, r = 0.773, p = 0.015, two-sided), together with the MFC loadings for each distance. Red line and shade represent the linear regression line with 95% confidence intervals. Source data are provided as a Source Data file. Permutation test is used by randomly shuffling the orders of regional MFC 5000 times to assess the association. Analyses and visualization were performed using custom scripts in MATLAB R2022b (MathWorks, Natick, MA, USA) together with Connectome Workbench.