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

Domain: chemistry
Published: 2026-03-19
Editorial outcome: major_revision

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.

Introduction
The perinatal period and early postnatal weeks constitute a critical epoch for cerebral development, characterized by orchestrated multidimensional reorganization of neural architecture and emergent functional networks1–4. During this phase, rapid cortical expansion is accompanied by explosive synaptogenesis, dendritic arborization, and myelination processes that collectively scaffold the nascent functional networks observed in neonates by 37-44 weeks of gestational age5–7. Since the first adult-like default-mode functional network was identified in neonates’ brain4, a hierarchical process of early brain functional network development has been delineated over the past two decades3,8. Functional networks develop through hierarchical circuit assembly: starting before birth, primary sensory/motor integration stabilizes by 3 months of age, limbic-cortical regulation emerges for initial emotional control by 24 months, and advanced associative networks like the default-mode network begin crystallizing after 2 years of age3,7. However, the underlying mechanisms linking discrete functional capacities with structural elements, in particular myelination, which establishes the biophysical substrate for signal propagation, remain elusive.

The coupling mechanisms9,10 of brain structure and function are key to understanding early brain functional integration. Current efforts to map structure-function coupling in the neonatal brain reveal fundamental paradoxes that challenge conventional neurodevelopmental frameworks. While adult studies demonstrate moderate spatial correspondence between structural connectivity and functional connectivity, neonatal data display that long-range connections often exhibit weaker structure-functional alignment (lower coupling)10–13. This characteristic is interpreted as a functional specification of the cerebral cortex during the early development12,14–16, allowing these regions to display more flexible functional connectivity that is not strictly constrained by physical axonal wiring. The distance-dependent decoupling mechanism contradict the hypothesis that functional networks emerge from progressive myelination of long-distance white matter tracts8. More notably, the age-dependent decoupling phenomenon of brain structure and function during early development10 also contradicts other studies reporting that the structural-functional coupling undergoes sustained enhancement from 30 weeks of gestational age through approximately 20 years of age9,14,17. This increased coupling is driven by the parallel growth of both integration and segregation between the two connectomes. However, an important neural structural basis, myelination, in the infant brain is often overlooked. Myelination is known to be activity-dependent, thus coordinated developments with functional connectivity are expected18–22. However, there is currently no previous study examining the direct, quantitative relationships between coordinated myelination and functional synchrony in the developing brain. Moreover, although white matter myelination has been more extensively studied, it remains unresolved whether it is the myelination of long-range tracts or the maturation of local gray matter that correlates more closely with the emergent functional connectivity. These represent critical knowledge gaps in better understanding infant structural-functional development.

As extensively reported, myelination is immature in the newborn stage. However, despite far-from-mature myelination levels in white matter tracts23,24 (30% of the adult level), neonates exhibit resting-state functional networks sharing 80% spatial similarity with adult counterparts25,26. This suggests that white matter-based structure-function coupling is likely insufficient to fully explain the underlying dynamics of the developing infant brain. Therefore, we hypothesize that the formation of functional connectivity does not solely rely on white matter fiber connections. Instead, we propose that coordinated growth of local gray matter myelination is more closely related to the emergence of initial functional connectivity. The brain wiring, including synaptic formation and pruning within the signaling-integrated white matter, is more potently modulated by the distributed gray matter systems27–31. Notably, gray matter, which primarily consists of neuronal cell bodies, develops earlier in the fetus than white matter, and undergoes rapid growth in late gestation and continues to develop after birth5–7. Unlike white matter myelination that peaks postnatally, gray matter oligodendrogenesis commences as early as 28 weeks of gestational age, synchronizing with the emergence of local field potential oscillations1,27,28,32,33. Such findings imply that early functional architectures may arise from the cumulative cortical growth.

Homophilic principles of nervous system organization state that brain regions with similar developmental trajectories will be similar also in terms of other properties33–35, hypothetically owing to shared genetic and/or convergent environmental effects. The cerebral cortex undergoes a process of spatially organized and heterochronic maturation34,36. Structural covariant connectivity characterizes the coordination of neural development measured by statistical similarity in developmental trajectories or distributions of structural indicators between paired brain regions. Over the past two decades, various methods of structural covariance analysis have been developed37–40, primarily encompassing both population-level and individual-level dimensions, and were used to study heterochronic network development41,42. Emerging evidence implicates gray matter myelination as a pivotal yet understudied mediator18,27,43,44. Cortical myelination is correlated with cortical development and intracortical circuit complexity45,46, and gray and white matter myelination are significantly correlated47. The first wave of myelination initiates prenatally, during which the primary sensory and motor cortices undergo myelination, and heteromodal association cortices do not become fully myelinated until the second decade, which is similar to the hierarchical development of cortical function5,18,48. Additionally, the structural covariance connectome also has similar dense characters49 with the functional connectome that even strong connectivity could be performed without direct fiber connectivity. Therefore, we held the hypothesis that the homophilic cortical growth measured by the myelination covariance could support the emergence of functional synchronization during infancy.

This study proposed a dual-axis myelination covariance framework to quantify the homophilic development, which encodes how regionally heterogeneous myelination schedules co-vary across both the cohort and neuroanatomical space, captured by two dimensions of covariance: (1) inter-subject covariance calculated as the group-level myelination covariance (gMC), representing synchronized myelination trajectories across subjects, and (2) intra-subject covariance calculated as subject-level myelination covariance (sMC), mapping coordinated myelination stages between distinct brain areas within individuals. This dual-axis covariance framework provides the statistical foundation for investigating spatially coordinated myelin growth that may show activity-dependent growth with emergent functional connectivity. Based on the dual-axis myelination covariance computed from myelin-sensitive MRI metric (T1w/T2w ratio), we further proposed a structure-function coupling approach called myelination-function coupling (MFC), which is specific for newborn brain development (Fig. 1) and is superior to traditional white matter-based structure-function analyses. A distinct MFC pattern was identified in neonates with multimodal brain imaging data from the Developing Human Connectome Project (dHCP)50. The single axis myelination covariance, either sMC or gMC, was also derived sMFC or gMFC to compare. Then, this study explored several factors associated with the MFC, including the age effects, the distance dependence, the birth effects, the transcriptomics bases, and the behavioral association. Results indicated that this coupling index can unveil the distinct features of how the brain's functional organization emerges from its parallel structural basis during the newborn stage.

Results
The dominant myelination-function coupling in the sensory cortex

We first examined the presence of association between functional connectivity (FC) and myelination covariance connectivity (MC) at both group and individual levels (Methods: Calculation of myelination-function coupling) using the public multi-modal imaging data for full-term neonates from the dHCP (Methods: Participants and data acquisition). The group level whole-brain connections of the two single-axis myelination covariance matrices were both significantly aligned with functional connectivity (gMC for the intra-hemisphere/whole connectome: R2=0.30/0.23, β=0.45/0.38, p<0.0001/p<0.0001, sMC for the intra-hemisphere/whole connectome: R2=0.22/0.23, β=0.097/0.087, p<0.0001/p<0.0001; Fig. 2a), which were consistent to the inter-hemisphere connections (see Supplementary Fig. 1a). The averaged maps of individual MFC, as well as those of gMFC and sMFC, were illustrated (Fig. 2b, c), showing an expanded higher MFC level than either gMFC or sMFC. Dominant regions in the MFC map included the primary motor cortex, posterior insular cortex, lateral temporal cortex, and primary visual cortex. Distinct regions—such as the primary auditory cortex and primary visual cortex—were dominant for gMFC, whereas regions including the primary motor cortex, insula cortex and inferior temporal cortex were mainly dominant for sMFC. Among the three coupling indices, gMFC and sMFC exhibited lower similarity (corr = 0.38, p<0.0001), which was substantially lower than their similarities with MFC (corr = 0.78 for gMFC and 0.83 for sMFC, both p<0.0001). Spin permutation tests on the three coupling indices revealed that the majority of vertices had significantly higher MFC than the null coupling level following 5000 randomizations (p <0.05; Supplementary Fig. 1b). Moreover, the coupling indices at the system level (Yeo’s 7-network parcellation51) also showed significance (Fig. 2d and Supplementary Table 1). The differences observed in mean values among the three coupling indices were further confirmed by independent two-sample t-tests at the functional system level (Supplementary Table 2).

To further investigate the spatial hierarchy of MFC, we computed the principal gradient of the FC connectome and correlated it with the MFC pattern. High consistency was observed (corr = 0.69 at Icosahedron-1002 atlas level, p<0.0001; Supplementary Fig. 1c; corr = 0.67 at vertex level, p<0.0001), suggesting that the MFC hierarchical pattern also captured the regions’ functional gradient variance. In addition, the proposed dual-axis MFC framework was repeated at the atlas level, yielding atlas-level MFC values using two alternative approaches for individual structural covariant connectome: KL divergence-based sMC (Kullback-Leibler divergence52,53 of probability distribution of myelin between paired brain regions) and the atlas-level sMC based on dot product proposed in this study. Using neonatal AAL atlas, comparisons among the vertex-level MFC map (after averaging to the atlas level) and the two atlas-level MFC maps revealed highly similar spatial patterns overall (corr = 0.52, p<0.0001 for KL divergence-based MFC, and corr = 0.61, p<0.0001 for dot-product-based atlas-level MFC, compared with the atlas-mean MFC; Supplementary Fig. 2a), with higher similarity observed between the two atlas-level MFCs (corr = 0.81, p<0.0001). Using Icosahedron-1002 atlas with fine-granularity, a more consistent MFC pattern with the vertex-level was further observed (corr = 0.98, p<0.0001 for vertex-level MFC and atlas-level MFC; Supplementary Fig. 2b). Considering the subcortical regions, the MFCs were highly consistent (corr = 0.98, p<0.001) between that either with or without the subcortical regions in the calculating framework (Supplementary Fig. 3 and Supplementary Result 1). Furthermore, a direct coupling between the functional connectivity profile and the myelin map showed significantly lower coupling than our proposed MFC (see Supplementary Fig. 4 and Supplementary Result 2).

Limited contribution from the white matter connectivity for the coupling

Having established that structural covariance of gray matter myelination is associated with the emergence of functional connectivity, it remained possible that white matter fiber connections also play a role therein. Therefore, we incorporated diffusion tensor imaging (DTI) data for fiber tract reconstruction, calculated the fiber tract connectivity density between brain regions to construct the structural connectivity network (SCN) based on the UNC infant AAL atlas. Subsequently, the fully weighted SCN was used to model the coupling relationship between white matter tracts and functional connectivity. The coupling was calculated using the same method as MFC, adopting a linear regression model. Averaged map of atlas-level individual MFC and SFC were illustrated (Supplementary Fig. 5a, b), showing that a moderate level of correlation between them (corr = 0.48, p<0.0001; Supplementary Fig. 5c). MFC map showed much higher level than SFC across the whole brain (the mean MFC is 0.22, and the mean SFC is 0.09), except for the cuneate cortex, with the peak increment up to 9 times (Supplementary Fig. 5d, the coupling increment map of MFC over SFC). The dominant regions of SFC were only distributed in the medial cortex, including the primary visual cortex and primary motor cortex. Furthermore, we introduced structural connectivity (SC) as an additional predictive variable into the MFC dual-axis framework. For the majority of brain regions, the introduction of SC profile barely improved the coupling strength, particularly in the dorsolateral prefrontal cortex (Supplementary Fig. 5e, MFC map with SC included). Regions with an increment exceeding 20% were concentrated in the posterior cingulate cortex, calcarine cortex, and cuneus (Supplementary Fig. 5f, Coupling increment map of SC-included MFC over MFC without SC variable).

The spatial growth patterns of myelination-function coupling

To identify the typical development trajectories of MFC during the neonatal period (Fig. 3a), we applied a generalized additive model (GAM) to capture both linear and non-linear relationships between MFC and postmenstrual age at scan (in weeks). Vertex-wise GAM for MFC demonstrated broadly distributed significant growth, including the whole visual cortex, inferior temporal cortex, parieto-occipital junction, sensorimotor cortex, insula, and prefrontal cortex (pFDR <0.05; Fig. 3b). Partially, the sensorimotor cortex, insula cortex, temporal cortex, and parieto-occipital junction exhibited significance for sMFC, whereas only the visual cortex was significant for gMFC (Fig. 3c). Notably, the primary motor cortex with dominant MFC at the vertex-level was not dominant in this growth pattern. Significant growth of whole-brain averaged MFC was also identified using GAM, with similar coefficients for sMFC and gMFC (β = 0.0049, p<0.0001 for both; Supplementary Fig. 6a), and a higher coefficient for MFC (β = 0.0083, p<0.0001). When parcellating the whole brain into seven systems and applying locally weighted scatterplot smoothing (LOWESS), the three coupling indices showed parallel increasing across all seven systems (Supplementary Fig. 6b), with the highest values in the sensorimotor system. Notably, sMFC and MFC showed a higher growth slope in the visual cortex than gMFC. The linear model was repeated for the global-level and vertex-level coupling indices to assess stable linear effects (Supplementary Fig. 6c, d). The growth-rate coefficients were delineated for each system, with the highest values in the visual and motor systems, and MFC exhibited the highest coefficients across all systems except the sensorimotor system, relative to sMFC and gMFC (Fig. 3d). Coefficient comparisons showed that sMFC was significantly higher than gMFC in the sensorimotor (SM), frontoparietal (FP), default mode (DMN), limbic (LIM) and ventral attention (VA) systems, while gMFC was higher in the visual (VIS) and dorsal attention (DA) systems (Fig. 3d and Supplementary Table 3). The spider plot (Fig. 3e) illustrated the percentage of significantly growing regions per system, indicating complementary distributions of sMFC and gMFC growth across the seven systems, and that MFC could capture all such growth effects.

The distance dependence of myelination-function coupling

Considering distance-related wiring costs during functional network emergence and potential overestimates of the coupling between myelination and functional profile caused by the spatially adjacent vertices, the cortical distance dependence of myelination-function coupling was further explored. The vertex-to-vertex connection profiles were ranked according to the geodesic distance between vertices and equally divided into multiple subsets from local to remote connections. Then, a series of MFC maps was calculated. The MFC spatial patterns exhibited high dynamics from the local to the remote connections, and each subset of connection showed distinct spatial MFC pattern (nine subsets, about 11% per range for each; Fig. 4a). Regions including the primary motor, primary visual and auditory cortex displayed dominant MFC for short-range connections, whereas the inferior temporal lobe and medial prefrontal cortex showed high MFC for long-range connections. Using finer range divisions (fifty subsets, 2% per range; the similarities between MFC maps of paired subsets in Supplementary Fig. 7a), the distance-dependent MFC curves within each system were further illustrated (solid curves in Fig. 4b), showing significantly higher MFC (around 0.4) for the most local connections and relatively lower values (around 0.15) for other distance ranges. Notably, the distance-dependent decrease in MFC exhibited an inflection point at around 12% percentile across all systems. Prior to this point, MFC decreased rapidly and linearly, whereas after this point, MFC varied smoothly and non-linearly from local to remote connections.

Leveraging these distance-dependent subsets, we further generated distance-dependent linear growth rate curves for each system (see dashed curves in Fig. 4b), showing a gradual increase from local connections (growth rate consistently around zero) to remote connections, with distinct patterns for each system. Linear relationships between the growth rate and geodesic distance were mapped as coefficients on the cortical surface (Fig. 4c). Positive regions indicated higher growth rate for long-range connections, as exemplified by distinct regions including the middle temporal gyrus, the premotor cortex and supplementary motor area (Fig. 4c). Negative regions indicated higher growth rate for short-range connections, as exemplified by the distinct regions including the primary motor, the posterior insula and the inferior parietal lobule (Fig. 4c). To evaluate the potential distance effect on the MFC patterns, we updated the MFC maps using two approaches: distance profile regressing (Fig. 4d) and local connections removing (Fig. 4e). Both maps showed high spatial similarity to the original MFC map, with significant correlations (corr = 0.86 and 0.88, p<0.0001; Supplementary Fig. 7b).

Similar distance-dependent MFC patterns were also reproduced across the five cyto-architectonic classes (Supplementary Fig. 7c) and for several critical regions (Supplementary Fig. 8). Convergently in details, the MFC of VIS and DA exhibited higher growth sensitivity (highest growth rate) to the medium distance connections (Fig. 4b), which was similar to that of the polar (POL) and parietal (PAR) cyto-architectonic classes; the MFC growth rate of DMN was lower for medium distance connections but higher for both local and remote connections; the MFC growth rate of VA showed higher for remote connections, which was consistent with the agranular (AGR) and granular (GRA) cyto-architectonic classes; the MFC growth rate of FP was higher for local connections, a pattern shared with the PAR cyto-architectonic class. To further validate our results, we replicated key findings using the updated MFC with distance regressed (see Supplementary Fig. 9). These replicated analyses included: (1) the global association between gMC/sMC and mFC; (2) the similarity between the updated MFCs and the original MFCs for 50 subsets; (3) the distance-dependent MFC curves and growth rate curves; (4) the global MFC growth rate and vertex-wise growth-rate coefficient of MFC with pFDR <0.05.

The specific effect of birth on myelination-function coupling

Next, we sought to explore the differential influences of gestational age (GA) at birth (defined as the weeks from fertilization to birth) and postnatal age (PNA) at scan (defined as the time from birth to MRI scan) on the development of cortical MFC. These two age metrics reflect the distinct effects of intrauterine and extrauterine environments, considering birth as a critical event involving environmental changes. We compared two prediction models for the vertex-wise MFC using different explanatory variables: (I) only postmenstrual age (PMA) at scan; (II) GA and PNA combined. The adjusted coefficient of determination R2adj of model II (based on GA and PNA combined) outperformed that of model I (based on PMA alone) (Fig. 5a), with vertices more widely distributed in the high R2adj value range. Specifically, the mean R2adj of model II was higher than that of model I: for model II, vertices with pFWE <0.05 account for 95.2%, with a mean R2adj of 0.0368 (range: 0.002–0.146) at these significant vertices; for model I, vertices with pFWE <0.05 account for 92.4%, with a mean R2adj of 0.0306 (range: 0.003–0.140) at these vertices. The model fitness of the two models across different systems was also compared (Fig. 5c and Supplementary Table 7; only vertices with fitting effect meeting pFWE <0.05 were included). Higher R2adj were observed for Model II (all p<10–5), especially in the SM, DMN, and FP systems.

To exemplify the effect of birth for MFC, we attempted to measure the different developmental patterns of MFC before and after birth, using the linear estimated coefficients of MFC for intrauterine and extrauterine development separately. A positive MFC growth associated with GA was observed across the whole cortex (Fig. 5b), particularly in the motor and visual cortices. For extrauterine development, all regions exhibited positive MFC growth except for the primary motor cortex and the posterior cingulate cortex (Fig. 5b). The coefficients comparison between GA and PNA suggested that the GA effect was significantly greater than PNA effect (all p<10-5; Fig. 5c and Supplementary Table 6), indicating a more profound influence from the intrauterine development. To further verify the birth effect on MFC development, we extracted two participant subsets: (1) full-term neonates with small PNA, representing intrauterine development, (2) preterm neonates with term-equivalent PMA and long PNA, representing extrauterine development (Fig. 5d). Both MFCs showed similar cortical patterns with higher values in the motor cortex. In comparison, full-term neonates exhibited significantly higher MFC levels (p<10-5; Fig. 5e, f), especially in the inferior parietal lobule (IPL), precuneus (PCUN), and superior frontal gyrus (SFG).

The transcriptomic association of myelination-function coupling

To explore the transcriptomic basis of the myelination-function coupling during neonates, this study utilized the BrainSpan Atlas of the Developing Human Brain, aged from 24 post-conceptual weeks (pcw) to 37 pcw with 11 regions. The first component of partial least squares (PLS) depicted gene expression pattern that was mostly associated with MFC distributions, which explained 20.38% of the variance in gene expression data (pperm = 0 with 5,000 times permutation tests; Fig. 6a). Meanwhile, under the PLS1 (first component of PLS), significant spatial association between the MFC score and the gene expression score was observed (corr = 0.91, p=0.0001; Fig. 6b). Furthermore, by utilizing a univariate one-sample Z test, we ranked these candidate genes based on normalized PLS1 weights, identifying 352 positively weighted genes (PLS1+: Z>3, pFDR <0.001; Fig. 6c) and 112 negatively (PLS1-: Z<−3, pFDR <0.001; Fig. 6c). We utilized the Metascape54 to align PLS1+ and PLS1- genes with the gene ontology (GO) terms, including biological process (BP), molecular function (MF) and cellular component (CC). The PLS1+ gene list exhibited significant enrichment in terms related to angiogenesis, transmembrane transport, neural regulation and immune response (the top 20 GO enriched clusters shown in Fig. 6d, orange bars, and Fig. 6e), such as “side of membrane”, “tube morphogenesis”, “transepithelial transport” and “gliogenesis”. The PLS1- genes were associated with neuronal morphogenesis and neural regulation (see the Fig. 6d, blue bars), such as “synaptic membrane”, “neuron projection cytoplasm” and “axon”.

The behavioral association of myelination-function coupling

Behavioral assessments from the Bayley III scale for dHCP subjects at approximately 18 months of age were associated with the individual MFC values in this study. Non-parametric permutation analysis (5000 times permutation) for null-hypothesis testing revealed a statistically significant population correspondence between the regional MFC and behavioral profiles. Based on PLS, whole-brain MFC and behavioral scores for PLS1 showed a significant correlation (adjusted R2=0.026, pperm = 0.006; Fig. 7a), with the negative behavior loadings and positive image loadings. Similarly, under the PLS1, whole-brain functional connectivity strength (FCS) and behavioral scores also showed significance (adjusted R2=0.047, pperm = 0.001; Fig. 7b), with corresponding behavior and imaging loadings. For functional system-level PLS analysis, remote-distance sensorimotor MFC was significantly associated with the motor behavioral domain (adjusted R2=0.045, pperm = 0.041; Fig. 7c). Notably, for the 9 distance subsets (see Fig. 4a), the distance-specific MFC-behavior relationships (adjusted R2 scores for PLS1) exhibited increasing association with the connection distance, with adjusted R2 values increasing from 0.020 for the local connections to 0.092 for the remote connections (Fig. 7d and Supplementary Fig. 10).

Methods
Participants and data acquisition

A total of 887 scan sessions from 783 neonatal subjects were collected from the Developing Human Connectome project (dHCP). Based on the datasets screened in Logan’s research75, we further excluded subjects with missing myelin maps, calculated with the T1w/T2w ratio. Finally, 447 scans (447 subjects) were included in the following analysis, in which 364 were term born (female/male = 167/197; GA: 39.65±1.32 weeks; PMA at scan: 40.80±1.96 weeks), and 83 were preterm-born (female/male = 34/49; GA: 32.08±3.47 weeks; PMA at scan: 40.73±2.09 weeks). The dHCP study was approved by the UK Health Research Authority (14/LO/1169), and written consent was obtained from the parents or legal guardians of the study subjects. All neuroimaging data were acquired at the Evelina Newborn Imaging Centre, Evelina London Children’s Hospital, using a 3-Tesla Philips Achieva system. Detailed information on data acquisition parameters is available on the dHCP website (http://www.developingconnectome.org/).

MRI data preprocessing

All anatomical and functional MRI data were minimally preprocessed with the details of the dHCP structural pipeline76, and the functional pipeline77. Briefly, the surface-based pipeline for structural images included superresolution reconstruction to obtain the 3D T1w/T2w volumes78, registration (from T1w to T2w), bias correction, brain extraction, segmentation79, surface extraction80, and surface registration81. Furthermore, for the anatomical data, T1w/T2w intensity values were sampled at the mid-thickness surface. We used the new version of the Multimodal Surface Matching (new MSM81,82) to perform a non-linear alignment of surfaces and T1w/T2w values of each neonate to the dHCP age-specific templates, which reside in the HCP-YA fs_LR 32k space and cover the age range from 28 to 40 weeks PMA75,83. Subsequently, data registered to age-specific templates were aligned to the 40-week PMA template, downsampled to a 5k_fs_LR mesh, and lightly smoothed with 1.5 mm full-width-half maximum (FWHM) using workbench84 to better simulate the spatial continuity.

For the rs-fMRI data, the minimally pre-processed 4D functional volume images, which were motion and distortion corrected (MCDC) with FSL EDDY, denoised based on spatial independent component analysis (sICA) and nuisance regression, temporally band-pass filtered (0.01−0.09 Hz) and spatial registration. Even though the potential head motion effects on fMRI signals had been carefully controlled in MCDC, sICA and motion parameters regression, we further computed the relationship between mean relative root mean square (RMS) displacement and functional connectivity (QC-FC correlations)85, resulting that the distribution of QC-FC correlations was zero-centered with no significance. The residual mean framewise displacements (mFDs) were further examined for age effects, and nothing significant was found. For the registration challenge since the fast development of infant brain, a two-step template alignment method was adopted, with a non-linear registration from native functional space to each age specific template80 (dHCP volume atlas, https://gin.g-node.org/BioMedIA/dhcp-volumetric-atlas-groupwise), which included templates for 36 to 44 weeks of postmenstrual age, and another linear alignment from each age specific template to the target 40-week PMA template (dHCP age specific template). The two transform files were estimated and applied to the functional volumes using Advanced Normalization Tools (ANTS, http://stnava.github.io/ANTs/). Finally, the cleaned rs-fMRI data were further projected onto the individual native cortical surface and then spatially registered and smoothed using surface pipelines.

Functional connectivity network

Functional connectivity (FC) between each pair of vertices was quantified as the Fisher-Z-transformed Pearson correlation between the pairwise BOLD time series. For each participant, an 8589×8589 (8589 is the total number of vertices in the left and right hemispheres after removing subcortical structures at the 5k surface resolution) weighted adjacency matrix encoded the functional connectivity network. In addition, another two FC matrices with different cortical granularities were generated for the Icosahedron-1002 atlas86,87 and the neonatal AAL atlas88. Removing the subcortical regions, the neonatal AAL atlas included 82 cortical regions, and the Icosahedron-1002 atlas included 923 uniform cortical regions within each hemisphere.

White matter-based structural connectivity network

Structural connectivity (SC) between each pair of regions was calculated as the fiber tract connectivity density using DTI imaging data for fiber tract reconstruction based on the ROI-level atlas. We utilized two alternative volume-based infant atlases: the UNC infant AAL atlas88, comprising 8 subcortical structures and 82 cortical structures. The reconstruction of the diffusion data was performed in native space with DSI-studio (https://dsi-studio.labsolver.org/). After the generalized q-sampling imaging (GQI) with a diffusion sampling length ratio of 1.25, the whole-brain fiber tracking was conducted with the parameters including the angular cutoff of 60 degrees, step size of 1.0 mm, minimum length of 30 mm, and maximum length of 300 mm. The whole-brain fiber tracking process was performed with the FACT algorithm until 1,000,000 streamlines were reconstructed for each individual. The fiber tract connectivity density was quantified as the number of streamlines normalized by volumes of paired regions, yielding a 90×90 (90 parcellations with 82 cortical regions and 8 subcortical regions) matrix encoded the structural connectivity network (SCN) for each participant. Then, Sparse structural connections were converted into a communicability matrix to avoid the loss of structural information, resulting in a fully weighted form of the SCN based on the communication model89,90.

Dual-axis myelination covariant network

In this study, we proposed a dual-axis myelination covariance framework to quantify the spatiotemporal homogeneity of myelin developmental dynamics during infancy. The vertex-wise myelin (T1w/T2w) was first Z-normalized to quantify myelination levels relative to whole-brain baselines within each individual. On this basis, the dual-axis MC was derived: (1) group-level myelination covariance (gMC) across subjects captures inter-subject covariance, which is defined as the Fisher-Z-transformed Pearson correlation between each pair of cortical vertices using Z-normalized myelin data connected as an N×L sequence across subjects (where N is the number of subjects, 364 for term-born and 83 for preterm-born; L is the number of vertices), yielding an 8589×8589 covariance matrix. (2) Subject-level myelination covariance (sMC) captures individual-level covariance patterns, defined as the intra-subject dot products of the Z-normalized myelin values, independent of confounding extraneous subject-level variables, yielding another 8589×8589 covariance matrix.

Calculation of myelination-function coupling

Vertex-wise connectivity profiles, represented as vectors of connectivity strength of a vertex to all other vertices in the whole brain, were extracted from each column of the myelination covariance connectome or functional connectome. Myelination-function coupling (MFC) was then measured using a linear regression model. For a given vertex, the predictive variable was the myelination covariance (MC) profile, and the functional connectivity (FC) profile was the dependent variable. As shown in the univariate model:

FC = b0 + b1 × MC, (1)

where the intercept b0 and regression coefficients b1 were estimated model parameters, MC was gMC profile or sMC profile corresponding to the calculation process of gMFC or sMFC. In the present study, we assumed that a linear model with a dual-axis myelination covariance matrix as multiple predictors captures more comprehensive and integrated information about myelination-function coupling. As shown in the bivariate model:

FC = b0 + b1 × gMC + b2 × sMC: (2)

The goodness of fit per vertex represent the vertex-level MFC, quantified as the adjusted coefficient of determination (R2adj). The R2adj from the univariate model was used in the quantification of gMFC or sMFC, and the R2adj from the bivariate model was used to quantify MFC. Vertex-wise indices of MFC were averaged across participants to create a mean vertex-level coupling map (Fig. 2c).

In general, the coefficient of determination (R2) tends to increase when the number of parameters included in the model increases. One strategy to account for this addition is to adjust the R2 measure based on the number of predictors. Specifically, we calculated R2adj for each brain region or vertex:

R2adj = 1 − (1 − R2) × (N − 1) / (N − p − 1), (3)

where N is the number of samples and p is the number of predictors in the model, and is equal to p = 2 for the dual-axis framework (bivariate model, including both sMC and gMC).

Age effects analysis of MFC

To further validate the development of cortical MFC in human infants, we considered age as a continuous variable using a generalized additive model (GAM)91 to flexibly investigate linear and nonlinear relationships between MFC and age. For each cortical vertex, the GAM model was defined as follows:

MFCi,subj = βi,0 + βi,1 × f(PMAsubj) + covs, (4)

where the MFC of participant subj was the dependent variable, i represented the vertex, f(:) was a smooth term and PMAsubj was the age of subj at scan. Importantly, thin plate regression splines were used for the smoothing basis, and the residual estimates of the maximum likelihood (REML) method was used to estimate nonlinearities, penalizing nonlinearity to avoid over-fitting the data92. The averaged first derivative of the age smooth function (ΔMFC/ΔPMA) characterized the age effects on MFC development. We controlled false positive errors for multiple comparisons using the False Discovery Rate (pFDR <0.05). The covs was co-variables including head size, sex and mean framewise displacement (mFD). In addition, linear age effects were also investigated using a linear regression model. Significance inference of linear effects was performed using surface-based, vertex-wise permutation testing (pFWE <0.05; FSL PALM93, version alpha 119) with threshold-free cluster enhancement (TFCE)94.

Distance dependence analysis of MFC

Considering that spatially adjacent vertices may intrinsically overestimate the coupling between myelination and functional profile, the three types of connection profiles (FC, gMC, sMC) were ranked according to the geodesic distance between vertices and equally separated into multiple subsets. This cortical geodesic distance between all pairwise vertices was calculated according to the workbench command “surface-geodesic-distance-all-to-all” on the mid-thickness surface, with values ranging from 0.75 mm to 139.11 mm (the mean distance for node subsets is detailed in Supplementary Table 8). Therefore, a series of MFCs was calculated to assess the coupling for the connections between nodes with various cortical distances, yielding a distance-specific examination for MFC. Detailed nine subsets (about 1/9=11% for each, mean value 18.47 mm for the most local) and fifty subsets (1/50=2% for each, mean value 8.12 mm for the most local) of the distances were selected to show distance-specific MFC patterns as well as the variation curves for systems or core regions. Furthermore, for each distance-specific MFC, the growth rates were assessed with PMA for each node to identify the growth sensitivity of distance-specific MFC. From the perspective of brain system integration, we further examined the distance-dependent properties of MFC within seven functional networks and five cytoarchitectonic systems. The MFC and its age effects for each system or network were derived by averaging across all constituent nodes.

Birth event effect analysis of MFC

To further explore the birth effects, we utilized a comparison between an early, experience-independent phase and a subsequent experience-dependent phase. Two distinct phases were separately represented by gestational age (GA) and postnatal age (PNA). PNA was the difference between GA at birth and the postmenstrual age (PMA) at scan. Therefore, two linear models were fitted at the vertex level for the MFC, with Model I for the PMA effect:

MFCi,subj = βi,0 + βi,1 × PMAsubj + covs, (5)

and Model II for GA and PNA effects, accounting for the development of extrauterine or intrauterine,

MFCi,subj = βi,0 + βi,1 × GAsubj + βi,1 × PNAsubj + covs, (6)

where the GAsubj and PNAsubj were the age at birth and the age at scan after birth. The birth effect was indicated by the difference between Model I and Model II, considering the adjusted R2 of each model and the linear coefficients of the age terms. To disentangle the influence of early postnatal experience from endogenous developmental processes, we compared the group mean MFC in two sub-cohorts matched on PMA: a preterm-born cohort (N=83; female/male = 34/49; GA: 32.08±3.47 weeks; PMA at scan: 40.73±2.09 weeks) and a matched term-born cohort (N=83; female/male = 38/45; GA: 39.65±1.32 weeks; PMA at scan: 40.80±1.96 weeks). As the preterm-born cohort had a shorter gestational age but longer postnatal experience at the time of scan, any differences in MFC patterns can be more confidently attributed to the effects of birth and ex-uterine experience. Significance inference of the linear modeling analysis was conducted with permutation testing (FSL PALM93) with TFCE94 under pFWE <0.05.

Spatial association with transcriptional gene expression

To elucidate the underlying biological mechanisms of myelination-function coupling during neonates, this study utilized transcriptomic data from donors aged 24 post-conceptual weeks (pcw) and 37 pcw obtained from the BrainSpan Atlas (https://www.brainspan.org/static/home) of the Developing Human Brain. Samples with insufficient regional coverage were excluded under the assumption that the BrainSpan cohort and our study populations share homogeneity. RNA-Seq expression averaged values were extracted from 11 brain regions, including the primary sensorimotor cortices (M1C, V1C, A1C, S1C) and other higher-order cortices (see the Supplementary Table 9 for details). Then, PLS (partial least squares) regression was employed to probe the potential relationship between MFC and transcriptional levels for all genes. PLS1 (first component of PLS) depicted the most correlated cortical expression map with respect to the MFC. Additionally, we performed a nonparametric permutation test by randomly shuffling the orders of regional MFC 5000 times to assess the statistical significance of PLS1. The bootstrapping technique was performed with 5000 bootstrap replicates to evaluate the contribution of each gene to the PLS1. The Z-score per gene was determined by dividing its weight by the corresponding bootstrap standard error, based on which all genes were ranked. Then, genes exhibiting Z>3 and pFDR <0.001 were assigned to the PLS1+ gene list, while those with Z<-3 and pFDR <0.001 were allocated to the PLS1- gene list. Subsequently, the PLS1+ and PLS1- gene lists were separately subjected to the Metascape (https://metascape.org/gp-/index.html) to identify enriched GO terms, including biological process (BP), molecular function (MF), and cellular component (CC).

Analysis of behavioral outcomes

To investigate associations between the proposed MFC and the neurodevelopmental outcomes, we analyzed the cortical measures from the neonatal data together with their follow-up behavior assessment at about 18 months of age. Vertex-wise MFC maps were parcellated into 923 uniform cortical regions for each hemisphere to generate the individual MFC feature vector, yielding a 1×923 MFC feature vector (averaging from the two hemispheres). To compare, another individual cortical measure, the functional connectivity strength (FCS), also extracted into a 1×923 FCS feature vector and linked with the behavioral outcomes. Neurodevelopmental outcomes were assessed using the Bayley Scales of Infant Development III, including 95,96 composite scores for: motor development (fine/gross motor skills), language development (receptive/expressive communication), and cognitive development (sensory perception, problem-solving). PLS regression95,96 was applied to the brain matrix (participants × 923 cortical measures) and behavior composite scores (participants × 3 Bayley domains). Significance of PLS1 was determined through 5000 times permutation test (pperm < 0.05) for multiple comparisons. To interpret the PLS1, we computed PLS loadings, measured by Pearson’s correlations between the original brain data and brain composite scores, as well as between the original behavior data and behavioral composite scores for PLS1. To estimate confidence intervals for these loadings, we applied a bootstrapping procedure that generated 500 samples from subjects’ brain and behavior data. Code for this analysis can be found here: https://github.com/MIPLabCH/myPLS.

Null model

We tested the observed coupling against the spatial permutation test (“spin test”). Using this null model, we explored whether the observed coupling was specific to the actual connectivity pattern rather than due to the spatial autocorrelation. Specifically, we first recorded the spherical coordinates of each vertex at a hemispherical resolution of 5k. Then, we randomly rotated the vertices while maintaining spatial autocorrelation and reassigned vertex values to the nearest vertices. This procedure was repeated 5000 times. We then calculated the MFC of null model to build a distribution of null maps that preserves spatial neighborhood information. The spin-based p-value was calculated as the proportion of MFC in the null model that exceeded (for positive coupling) or were weaker than (for negative coupling) the empirical MFC.

Connectome gradients

To identify spatial axes in connectivity variation across different areas, we applied an unsupervised technique called diffusion map embedding to compute gradients of the three connectomes separately (mFC, gMC, sMC). This technique converted the connectivity matrix into a normalized angle matrix, which scales the angle between each pair of vertices as a function of similarity. Manifold learning parameters were identical to those previously described, specifically α=0.05 and automated diffusion time estimation97. The computation was done with the MATLAB BrainSpace package (https://www.mathworks.com/matlabcentral/fileexchange/86887-brainspace).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This work investigates the relationship between cortical myelinations and functional connectivity in infants and proposes the concept that dual-axis myelination covariance drives the emergence of functional connectivity during infancy. MRI datasets from the developing human connectome project were used (including a large cohort of healthy neonatal subjects), which were processed according to standard pipelines. The authors mapped dominant myelination-function coupling (MFC) in sensor cortex and spatial growth patterns of MFC, examined cortical distance dependences of MFC and specificity effects of birth; transcriptomic and behavioral associations of MFC were also briefly studied. Overall the paper is well written, and the study is reasonably comprehensive. The reviewer has the following comments:\n1. In the section “The cortical distance dependence of myelination-function coupling”, Many regions were reported to have distance dependence of MFC. Readers may wonder if the distance dependence effect is a confound that should be controlled in MFC mapping.\n2. In lines 305-307: The adjusted R-squared for both model I and II appear to be rather low with no p-value reported, which raises the question about the reliability of the model comparison.\n3. In lines 824-842: The construction of dual-axis myelination covariant network is not quite clearly presented. It is recommended to expand Figure 1 with greater details so that readers can better understand the matrix construction process.", "reviewer_2": "In this work, Liu et al. explore links between functional connectivity and gray matter myelination using the dHCP dataset. The authors report that myelin development and the development of functional connectivity follow overlapping spatial patterns. I find this an interesting observation. At the same time, I had some difficulty following parts of the manuscript, and I believe some of the interpretations are not supported by the data. Below, I outline my main concerns:\n1. Throughout the manuscript, the authors conclude that local myelination, rather than white matter connections, may support functional connectivity during development. I struggle with this interpretation. By definition, functional connectivity between regions requires structural connections, so it is difficult to reconcile the claim that local myelination could support connectivity in isolation. Given emerging evidence that gray and white matter myelination are linked, I suspect the authors’ findings may instead reflect related changes in underlying white matter tracts.\n2. Since myelination is known to be activity-dependent, coordinated development of myelin and functional activity is expected. I find it problematic that the authors interpret myelination changes themselves as directly supporting functional synchrony.\n3. The interpretation of the T1w/T2w ratio as a measure of myelination is still debated. While more specific measures (e.g., R1 or MWF) are not available in the dHCP, I think this limitation needs to be explicitly acknowledged and discussed in more depth.\n4. To my understanding, the authors assess covariance in functional activity and covariance in myelination. I find it difficult to follow the interpretation of this as “connectivity.” Could the authors clarify what is meant by “myelination connectivity”?\n5. The reported association between functional connectivity and myelination seems to be driven largely by neighboring vertices. Could their main observations simply be a consequence of the fact that both measures change smoothly in space?\n6. In the transcriptomic analyses, I was surprised that no genes related to myelination were identified. This absence seems important and deserves discussion.\n7. For the behavioral associations, it would be interesting to examine whether networks predict their expected functional domains (e.g., whether motor system MFC predicts motor ability).\n8. Please include p-values in the main text.\n9. Could the authors explain why T1w/T2w measures were z-scored? I am concerned that this procedure might obscure meaningful inter-subject differences, such as those related to age.\nI did not try to run the code, I just looked at it briefly. It looks well organized and easy to follow. However, it seems not all anaylses are included at the moment? E.g. I could not find code to run the behavioral associations.", "reviewer_3": "This paper investigates the relationship between myelination and functional connectivity (FC) during infancy using fMRI data. While the topic is of interest, several significant methodological and analytical concerns need to be addressed.\nMajor Concerns:\n1. Head Movement Artifacts in Infant fMRI: Infants are prone to significant movement, and the manuscript does not provide a robust description of the head movement correction pipeline. Large head movements can introduce spurious functional correlations that could confound the central findings.\n2. Exclusion of Subcortical Structures: Early brain development is characterized by the maturation of deep-seated subcortical structures and their connections to the cortex. For example, the early development of connections between the thalamus and cortex might be a covariant cause for both functional connectivity and myelination development. The observed covariance between myelination and FC could be driven by these subcortical connections, and their exclusion from the analysis could lead to an incomplete and potentially misleading conclusion.\n3. Novelty and Differentiation is limited: The reported patterns of \"coupled regions\" in the results appear to be broadly consistent with previous literature that sensory systems are developed earlier. It is difficult to conclude that the findings are uniquely due to myelination or other well-known developmental factors. The analysis would be significantly strengthened by leveraging the Diffusion Tensor Imaging (DTI) data available in the Developing Human Connectome Project (dHCP) dataset to also model white matter tract development and its relationship with the observed functional changes.\n4. The relationship between MFC and Distance: The manuscript notes a strong presence of MFC in short-distance connections. The paper does not adequately explore whether this is due to a lack of myelination in white matter, which would imply that myelination is more critical for long-range FC development. It is highly suggested that distance be defined based on fiber length or myelination level, rather than a simple geodesic distance, to provide more meaningful insights into this relationship.\n5. Weak Behavioral Association: The effect size of the behavioral association is reported to be too low, which significantly weakens the clinical and biological relevance of the findings. The manuscript needs to provide a more compelling argument for the significance of this association despite the low effect size.\nMinor Concerns:\n1. R-squared Correction: It is not clear from the manuscript if any correction was done for the R2 results, which is a standard procedure in such analyses to account for the number of predictors."}, "author_response": "Dear Editor and Reviewers,\nWe sincerely thank the editor and the reviewers for their time and insightful feedback, which has been invaluable in strengthening the rigor, clarity, and impact of our study. During this revision, we have carefully considered the constructive comments raised by the Reviewers. And accordingly, we have thoroughly and thoughtfully addressed every point raised in this response letter.\nThe major revisions and enhancements include:\n1. Direct integration of white matter analysis to clarify gray matter contribution.\n2. Deepened analysis and discussion of distance effects.\n3. Expanded physiological interpretation of gene findings and discussion.\n4. Strengthened specificity and mechanistic interpretation of behavioral associations.\n5. Enhanced overall clarity, rigor and gramma.\nFor convenience in checking the changes made in the revised manuscript, we have marked these changes with ‘red’ color. We wish our revised manuscript would be more readable than the previous version and reach the high standard requirement of this journal.\nThank you once again for your time and expertise, and we are looking forward to your insightful comments again.\nBest wishes,\nYours sincerely,\nWeijin Liu & Yuanyuan Chen\nOn behalf of all authors\nREVIEWER COMMENTS AND AUTHOR RESPONSE\nReviewer #1 (Remarks to the Author):\nThis work investigates the relationship between cortical myelinations and functional connectivity in infants and proposes the concept that dual-axis myelination covariance drives the emergence of functional connectivity during infancy. MRI datasets from the developing human connectome project were used (including a large cohort of healthy neonatal subjects), which were processed according to standard pipelines. The authors mapped dominant myelination-function coupling (MFC) in sensor cortex and spatial growth patterns of MFC, examined cortical distance dependences of MFC and specificity effects of birth; transcriptomic and behavioral associations of MFC were also briefly studied. Over all the paper is well written, and the study is reasonably comprehensive. The reviewer has the following comments.\nResponse: Thanks for your efforts and kindly reviewing our manuscript. During these days, we have carefully considered and responded to the following concerns. Please kindly refer to these replies for modifications made in our revised manuscript.\nComment 1. In the section “The cortical distance dependence of myelination-function coupling”, many regions were reported to have a distance dependence of MFC. Readers may wonder if the distance dependence effect is a confound that should be controlled in MFC mapping.\nResponse: We thank you for raising this important methodological concern regarding distance as a potential confounding effect in our proposed MFC mapping. We fully agree that regressing out distance is crucial to ensure that the observed coupling reflects genuine biological coordination beyond mere spatial autocorrelation.\nIn our study, we have already addressed this concern. Specifically, following the established practices in the field,1,2 we regressed out the shared variance associated with geodesic distance from both the myelination covariance and functional connectivity profiles before computing the coupling index (MFC) at each vertex. Subsequently, the spatial correlation between MFC maps with and without regressing out the geodesic distance was computed, as presented in the manuscript (Figure 4d and S. Figure 7b). These two maps demonstrate strong spatial consistency (corr=0.86, p< 0.0001), thereby validating that the global spatial pattern of our primary finding - the alignment between myelination covariance and functional connectivity - is robust and not driven artifactually by spatial distance. Here, we provided the revised Figure 4d and the Supplementary Figure 7b. For the revision of Figure 4b, we regressed out the distance for both functional connectivity and myelination covariance profiles, while only the functional connectivity profile was controlled for distance in the previous version.\nFigure 4d. Averaged map of individual MFC with distance regressed.\nSupplementary Figure 7b. The relationship between the original MFC and the MFC with regressing out the distance.\nIn addition, to visually present the distance-dependent pattern in the original manuscript, we have performed the distance dependence analysis for each vertex FC. The connections to all other vertices were equally divided into fifty subsets from local to remote, and a series of the distance-specific MFCs was obtained (as presented in Figure 4b, and Figure 4a for nine subsets version). Although these results confirmed the high coupling for the local connections, we further indicated the specific distance dependence for each region or network, like SM and LM in Figure 4b. Within each distance subset, most of the MFC is not zero without the global spatial consequence.\nYou astutely note the reported distance dependence within local connection subsets. We wish to clarify that this is a distinct and complementary finding. The strong global pattern persists after distance correction (corr=0.86, p< 0.0001; Supplementary Fig. 7b), while the analysis of connection subsets reveals how the strength of MFC itself varies systematically with connection length (Figure 4a illustrates the distinct spatial pattern of MFC for each distance subset). These results indicated a biologically informative characteristic of local network organization that does not invalidate the primary, distance-corrected result.\nTo further directly address your concern, we have performed three additional validation analyses during revision.\nFirst, we re-analyzed the global coupling between both MCs and FC: At the intra-hemisphere level, the two single-axis myelination covariant matrices (gMC and sMC) both remained significantly correlated with the functional connectivity, as shown in Figure R1a.\nSecond, the persistence of distance dependence pattern was confirmed: as displayed in Figure R1b, MFC calculated at fifty subsets of connections with various geodesic distances between one and all other vertices was unchanged, indicated by high correlation (corr=0.99, p<0.0001) between that of the original vertex-level MFC (8589 vertices × 50 subsets) and adjusted MFC (8589 vertices × 50 subsets) with controlling for distance. Therefore, it was reasonable that the distance dependence remained unchanged after regressing the distance (see Figure R1c, the distance dependence pattern for each network). This confirms that the described distance dependence is a stable feature of the data, not an artefact of the distance confounding effect that we already globally controlled for.\nFigure R1. Distance-dependent pattern of MFC after controlling for distance. a Linear regression between mFC and gMC or sMC for the connections at the connectome level. b Pearson correlation between the original vertex-level MFC of fifty subsets (8589 vertices per subset) and adjusted MFC with regressing out the distance of fifty subsets (8589 vertices per subset) is extremely high (r=0.99, p<0.0001). c After controlling for distance, 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 characteristics within the right and left hemispheres.\nThird, the developmental trajectories were examined after distance controlling: we recalculated the MFC growth rate derived from a generalized additive model (GAM) after controlling for distance. The significant developmental growth of MFC with age remains virtually identical during the neonatal period after regressing out distance (Figure R2a). The spatial map of significant growth regions shows high overlap with the original results (Dice coefficient = 0.79; see Figure R2b for results with controlling distance and see Figure 3b for original results).\nFigure R2. The spatial growth pattern of MFC after controlling for distance. a After controlling for distance, GAM fitting of MFC global values and PMA (p<0.0001). b After controlling for distance, growth-rate coefficient of MFC during infancy derived from GAM analysis with pFDR<0.05.\nIn summary, we do agree that the spatial distance is crucial to setup the functional path, and we still choose to include the MFC without distance correction in the subsequent analysis, primarily based on the three findings: first, with or without the distance control, the pattern of MFC did no change indicating its little impact on the MFC; second, the MFC pattern showed distinct distance dependent patterns; third, the MFC for longer connections presented higher value and association with behavior. The spatial distance is a double-edged confound, and it provides distinct biological characteristics about how the early brain function organizes with our distance dependence analysis.\nWe hope this revision addresses your concerns and clarifies the scientific question and motivation of our current work. If more modifications are still needed, please contact us to resolve this again.\nAnd we have clarified these points as follows in the revised manuscript (section “The distance dependence of myelination-function coupling” of the Results, line 284 – line 288 in manuscript):\n“To evaluate the potential distance effect on the MFC patterns, we updated the MFC maps using two approaches: distance profile regressing (Fig.4d) and local connections removing (Fig.4e). Both maps showed high spatial similarity to the original MFC map, with significant correlations (corr=0.86 and 0.88, p<0.0001; Supplementary Fig. 7b).\n… To further validate our results, we replicated key findings using the updated MFC with distance regressed (see Supplementary Fig. 9). These replicated analyses included: (1) the global association between gMC/sMC and mFC; (2) the similarity between the updated MFCs and the original MFCs for 50 subsets; (3) the distance-dependent MFC curves and growth rate curves; (4) the global MFC growth rate and its significant regions.”\nComment 2. In lines 305-307: The adjusted R-squared for both model I and model II appear to be rather low with no p-value reported, which raises the question about the reliability of the model comparison.\nResponse: We thank you for highlighting the importance of statistical rigor. We acknowledge that the mean R2 values of both models are low. This is because adj individual cortical locations demonstrate subtle and heterogeneous effects, and the reported mean value aggregates these effects, including many vertices with minimal relationship with age.\nCritically, R2 values of the vast majority of vertices are highly statistically significant. As we added in the revised manuscript, after rigorous family-wise error (FWE) correction for multiple comparisons across vertices, 95.2% of cortical vertices showed a significant fitting for Model II (based on GA and PNA combined), and 92.4% for Model I (based on PMA alone). This widespread significance confirms that the observed relationships are reliable and not due to chance.\nUnder this context of high statistical reliability, direct comparison between the two models is valid and informative to assess whether separating the effects of GA and PNA provides a better fit than using PMA alone. The consistently higher R2 values and greater extent of significant vertices for Model II robustly adj support our conclusion regarding the differential influence of intrauterine and extrauterine development on MFC.\nWe have also added the P-value and the adjusted R-squared range (Results: section “The specific effect of birth on myelination-function coupling”) in the revised main text (line 317 – line 324 in manuscript):\n“Specifically, the mean R2 of model II was higher than that of model I: for adj model II, vertices with p <0.05 account for 95.2%, with a mean R2 of 0.0368 FWE adj (range: 0.002 - 0.146) at these significant vertices; for model I, vertices with p <0.05 account for 92.4%, with a mean R2 of 0.0306 (range: 0.003 - 0.140) FWE adj at these vertices. The model fitness of the two models across different systems was also compared (Fig. 5c and Supplementary Table 4; only vertices with p <0.05 were included). Higher R2 were observed for Model II (all p<10-5), FWE adj especially in the SM, DMN, and FP systems.”\nComment 3. In lines 824-842: The construction of the dual-axis myelination covariant network is not quite clearly presented. It is recommended to expand Figure 1 with greater details so that readers can better understand the matrix construction process.\nResponse: Sorry for the unclear illustration and statements for Figure 1 and the main pipeline of the proposed work. We have expanded Figure 1 as follows, including detailed steps of the proposed myelination-function coupling (MFC) based on dual-axis covariance. As presented in Figure 1, there are four sections, including the question, brain data, MFC calculation, and main analysis of the MFC maps.\nSpecifically, for the calculation of MFC, the individual functional connectivity profile is explained by two myelination covariance (MC) profiles: one for the intra-axis (intra-subject axis) and the other for the inter-axis (inter-subject axis). The intra-axis is represented by the subject-level myelination covariance/similarity network (sMC)， calculated as the outer product of the z-scored cortical myelination vector. We also provided the alternative result with KL divergence-based sMC, in the Supplementary Fig. 2. The inter-axis is represented by the group-level myelination covariance (gMC), which is a typical structural similarity network. Therefore, when calculating the individual MFC map, the individual sMC and global gMC are used to predict the individual FC using the same columns. The fit goodness quantified as the adjusted coefficient of determination (R2 ) identifies MFC at each vertex, yielding the individual adj MFC map on the surface.\nWe have updated Figure 1 as: “\nFigure 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.”\nReviewer #2 (Remarks to the Author):\nIn this work, Liu et al. explore links between functional connectivity and gray matter myelination using the dHCP dataset. The authors report that myelin development and the development of functional connectivity follow overlapping spatial patterns. I find this an interesting observation. At the same time, I had some difficulty following parts of the manuscript, and I believe some of the interpretations are not supported by the data. Below, I outline my main concerns.\nResponse: Thanks for reviewing and providing constructive suggestions to improve our work. In this reply letter, as well as our revised manuscript, we have carefully replied to the following comments and tried our best to revise our manuscript. Please refer to the following one-by-one replies for more details.\nComment 1. Throughout the manuscript, the authors conclude that local myelination, rather than white matter connections, may support functional connectivity during development. I struggle with this interpretation. By definition, functional connectivity between regions requires structural connections, so it is difficult to reconcile the claim that local myelination could support connectivity in isolation. Given emerging evidence that gray and white matter myelination are linked, I suspect the authors’ findings may instead reflect related changes in underlying white matter tracts.\nResponse: We thank you for raising this insightful and fundamental comment, which allows us to clarify the conceptual framework and intended contribution of our study. We fully agree that the functional connectivity ultimately requires structural connections, and we do not argue that local myelination operates in isolation to support long-range functional connectivity. Actually, our core finding is complementary, not alternative. We seek to investigate whether, in the context of the profound immaturity of long-range white matter tracts in the neonatal brain, the coordinated maturational state of cortical regions (as myelination covariance in this study) provides a parallel and influential factor for the emergence and spatial patterning of functional networks during infancy.\nOur hypothesis is based on the developmental facts of gray matter function and white matter fibers in the neonatal period and the homophilic principles of nervous system organization:\n(1) The developmental asynchrony between gray matter and white matter in the neonatal period: As noted in our manuscript, myelination is immature in the newborn stage. However, despite far-from-mature myelination levels in white matter tracts3,4 (30% of the adult level), neonates exhibit resting-state functional networks sharing 80% spatial similarity with adult counterparts5,6. This suggests that white matter-based structure-function coupling is likely insufficient to fully explain the underlying dynamics of the developing infant brain (at least it appears so in the test data, ignoring technical flaws for the neonatal brain). Therefore, we hypothesize that the formation of functional connectivity does not solely rely on white matter fiber connections.\n(2) The homophilic principle of brain development: Homophilic principles state that brain regions with similar developmental trajectories will be similar also in terms of other properties7,8, hypothetically owing to shared genetic and/or convergent environmental effects. The feasibility of using structural covariance analysis to study heterochronic network development is anticipated by the finding that structurally similar and/or covarying brain regions tend to follow similar developmental trajectories9, and more synchronized neural activities.\nGiven that cortical myelination is correlated with cortical development and intracortical circuit complexity10,11, and that gray and white matter myelination are significantly correlated12, we held the hypothesis that the homophilic cortical growth measured by the myelination covariance during development could support the emergence of functional synchronization during infancy.\nTo make this assumption clearer for readers, we have rewritten the Introduction (line 79 – line 132 in manuscript):\n“However, an important neural structural basis, myelination, in the infant brain is often overlooked. Myelination is known to be activity-dependent, thus coordinated developments with functional connectivity are expected 21–23. However, there is currently no previous study examining the direct, quantitative relationships between coordinated myelination and functional synchrony in the developing brain. Moreover, although white matter myelination has been more extensively studied, it remains unresolved whether it is the myelination of long-range tracts or the maturation of local gray matter that coordinates more closely with the emergent functional connectivity. These represent critical knowledge gaps in better understanding of infant structural-functional development.\nAs extensively reported, myelination is immature in the newborn stage. However, despite far-from-mature myelination levels in white matter tracts 24,25 (30% of the adult level), neonates exhibit resting-state functional networks sharing 80% spatial similarity with adult counterparts26,27. This suggests that white matter-based structure-function coupling is likely insufficient to fully explain the underlying dynamics of the developing infant brain. Therefore, we hypothesize that the formation of functional connectivity does not solely rely on white matter fiber connections. Instead, we propose that coordinated growth of local gray matter myelination is more closely related to the emergence of initial functional connectivity. The brain wiring, including synaptic formation and pruning within the signaling-integrated white matter, is more potently modulated by the distributed gray matter systems 28–32. Notably, gray matter, which primarily consists of neuronal cell bodies, develops earlier in the fetus than white matter, and undergoes rapid growth in late gestation and continues to develop after birth.5–7 Unlike white matter myelination that peaks postnatally, gray matter oligodendrogenesis commences as early as 28 weeks of gestational age, synchronizing with the emergence of local field potential oscillations1,28,29,33,34. Such findings imply that early functional architectures may arise from the cumulative cortical growth.\nHomophilic principles of nervous system organization state that brain regions with similar developmental trajectories will be similar also in terms of other properties34–36, hypothetically owing to shared genetic and/or convergent environmental effects. Based on the known neurodevelopmental process of spatially organized and heterochronic maturation of the human cortex35,37, structural covariant connectivity characterizes the coordination of neural development measured by statistical similarity in developmental trajectories or distributions of structural indicators between paired brain regions. Over the past two decades, various methods of structural covariance analysis have been developed 38–41, primarily encompassing both population-level and individual-level dimensions, and were used to study heterochronic network development42,43. Emerging evidence implicates gray matter myelination as a pivotal yet understudied mediator21,28,44,45. Cortical myelination is correlated with cortical development and intracortical circuit complexity46,47, and gray and white matter myelination are significantly correlated48. The first wave of myelination initiates prenatally, during which the primary sensory and motor cortices undergo myelination, and heteromodal association cortices do not become fully myelinated until the second decade, which is similar to the hierarchical development of cortical function 5,21,49. Additionally, the structural covariance connectome also has similar dense characters 50 with the functional connectome that even strong connectivity could be performed without direct fiber connectivity. Therefore, we held the hypothesis that the homophilic cortical growth measured by the myelination covariance during development could support the emergence of functional synchronization during infancy.”\nMore importantly, to directly address your concern about the impact from white matter tracts, we further incorporated the analysis of white matter connectivity, which is also consistent with Comment #3 of Reviewer #3. To investigate whether our proposed MFC reflects changes in white matter tracts as you mentioned, we reconstructed white matter tracts to create a structural connectivity network (SC) and calculated a parallel white matter-based structure-function coupling (SFC) index using a similar methodology as that for our myelination function coupling (MFC). The key results are as follows:\nFigure R3. Comparison between spatial patterns of atlas-level SFC and MFC. a Averaged map of individual atlas-level SFC based on the neonate AAL atlas. b Averaged map of individual atlas-level MFC based on the neonate AAL atlas. c Spatial association between atlas-level SFC and atlas-level MFC (r=0.48, p<0.0001). d Coupling increment map of atlas-level MFC over SFC. e MFC with SC included (SC as an additional predictive variable). f Coupling increment map corresponding to MFC with SC included over MFC without SC variable shown in panel b.\nFirst, higher global explanation (coupling) was found in MFC than in SFC. We incorporated DTI data for fiber tract reconstruction, calculated the fiber tract connectivity density between brain regions (the number of streamlines normalized by the volumes of paired regions), and constructed the structural connectivity network (SCN) based on the ROI-level atlas (UNC infant AAL atlas, optimized for neonatal brain). Sparse structural connections were converted into a communicability matrix to avoid the loss of structural information, resulting in a fully weighted form of the SCN based on the communication model13,14. Subsequently, the fully weighted SCN was used to model the coupling relationship between white matter tracts and functional connectivity. The coupling was calculated using the same method as MFC, adopting a linear regression model. An averaged map of individual SFC at atlas level was shown as Figure R3a. The dominant regions of SFC were only distributed in the medial cortex, including the primary visual cortex and primary motor cortex. To conduct comparisons at the same level, we also calculated the MFC at the atlas level based on this infant AAL atlas (Figure R3b). A moderate level of correlation was observed between them (corr=0.48, P-value<0.0001, Figure R3c). MFC map showed much higher level than SFC across the whole brain (the mean MFC is 0.22, and the mean SFC is 0.09), except for the cuneate cortex, with the peak increment up to 9 times (Figure R3d, the coupling increment map of MFC over SFC), supporting our hypothesis that gray matter myelination covariance is more relevant to emerging functional connectivity during infancy.\nSecond, a little additive contribution was detected when including SC in the MFC framework. We added structural connectivity (SC) as an additional predictive variable into the MFC dual-axis framework. For the majority of brain regions, the introduction of SC profile barely improved the prediction strength for functional connectivity, particularly in the dorsolateral prefrontal cortex (Figure R3e, MFC map with SC included). Regions with an increment exceeding 20% were concentrated in the posterior cingulate cortex, calcarine cortex, and cuneus (Figure R3f, Coupling increment map of SC-included MFC over MFC without SC variable). Consistent results were observed between the alternative INFANT atlas15 and the neonatal AAL atlas (Figure R4).\nFigure R4. The spatial patterns of atlas-level SFC (IFANT atlas) and MFC. a. Spatial association between SFC and MFC with the INFANT atlas (r=0.70, p<0.0001). b Spatial association between atlas-level MFC and MFC with SC included (r=0.95, p<0.0001). c. Coupling increment percentage distribution of MFC with SC included over MFC without SC.\nAll these findings suggest that both white and gray matter myelination are correlated with functional connectivity, but the covariance of gray matter myelination across regions during infancy is more closely related to the emergence of functional connectivity. The observed myelination-function coupling (MFC) is not merely a proxy for contemporaneous white matter tract strength but rather captures unique variance in functional network organization. White matter-based SFC captures similar patterns, but with much lower levels of association, and when adding the white matter-based SFC into our framework, little additional information is found. This might indicate that our MFC pattern contains the information from white matter-based SFC. The lower degree of association between functional connectivity and major white matter myelination is actually expected since each major white matter tract is likely shared by many functional connections, thus one-to-one association is diminished. Actually, as you mentioned, the myelination of gray and white matter is significantly correlated.12 However, the principles and nature of white matter myelination of the human brain during early infancy are not well understood. One of the main theories about white matter myelin development during infancy proposes that postnatal myelination progresses in a spatially organized manner16,17. Different spatial gradients of myelination have been proposed, including that white matter myelination begins closest to the neurons and follows the direction of information flow4,18, or that it occurs along a proximal to distal axis across the brain. From this perspective, considering developmental immaturity in white matter, the local gray matter dominance during infancy also makes “economic” sense in that cortical functional coupling initially coordinates more with the less resource-demanding intracortical myelination before different regional processes propagate to the coordinated myelination of long-range major white matter pathways that need considerably more resources, time, and coordination.\nWe are grateful for your comment, which has pushed us to strengthen our evidence and clarify our model, significantly enhancing the conceptual contribution of our work. We hope this revision addresses your concerns and clarifies the scientific question and motivation of our current work. If more modifications are still needed, please contact us to resolve this again.\nTo make it clear for readers, we added details in the Results (section “The dominant myelination function coupling in the sensory cortex”, line 206 – line 211 in manuscript) and Supplementary Result 1 (see the Supplementary Information):\nIn the Results, “Regarding the white matter effect, we calculated the white matter-based structure-function coupling (SFC) map and it showed a similar spatial pattern (corr=0.48, p<0.001) with our proposed MFC map, but about eight times lower level than MFC (see Supplementary Fig. 3 and Supplementary Result 1). Including structure connectivity (SC) profiles in our MFC framework only gains limited benefits to the coupling level.”\nIn the Supplementary Result 1, “Having established that structural covariance of gray matter myelination is associated with the emergence of functional connectivity, it remained possible that white matter fiber connections also play a role therein. Therefore, we incorporated DTI data for fiber tract reconstruction, calculated the fiber tract connectivity density between brain regions to construct the structural connectivity network (SCN) based on the UNC infant AAL atlas (see Supplementary Method 1 for details). Subsequently, the fully weighted SCN was used to model the coupling relationship between white matter tracts and functional connectivity. The coupling was calculated using the same method as MFC, adopting a linear regression model. Averaged map of atlas-level individual MFC and SFC were illustrated (Supplementary Figure 3a and b), showing that a moderate level of correlation between them (corr=0.48, p<0.0001; Supplementary Fig. 3c). MFC map showed much higher level than SFC across the whole brain (the mean MFC is 0.22, and the mean SFC is 0.09), except for the cuneate cortex, with the peak increment up to 9 times (Supplementary Fig. 3d, the coupling increment map of MFC over SFC). The dominant regions of SFC were only distributed in the medial cortex, including the primary visual cortex and primary motor cortex. Furthermore, we introduced structural connectivity (SC) as an additional predictive variable into the MFC dual-axis framework. For the majority of brain regions, the introduction of SC profile barely improved the prediction strength for functional connectivity, particularly in the dorsolateral prefrontal cortex (Supplementary Fig. 3e, MFC map with SC included). Regions with an increment exceeding 20% were concentrated in the posterior cingulate cortex, calcarine cortex, and cuneus (Supplementary Fig. 3f, Coupling increment map of SC-included MFC over MFC without SC variable).”\nComment 2. Since myelination is known to be activity-dependent, coordinated development of myelin and functional activity is expected. I find it problematic that the authors interpret myelination changes themselves as directly supporting functional synchrony.\nResponse: Sorry for our unclear and confusing description. We totally agree with you that myelination is activity-dependent, and coordinated development of myelin and functional connectivity is thus expected. However, there is currently no previous study examining the direct, quantitative relationships between coordinated myelination and functional synchrony in the developing brain. However, we seek to investigate whether, in the context of the profound immaturity of long-range white matter tracts in the neonatal brain, the coordinated maturational state of cortical regions (as indexed by myelination covariance in this study) provides a parallel and influential factor for the emergence and spatial patterning of functional networks during infancy. As methodologically presented in the updated Figure 1 with the full pipeline of our proposed MFC calculation, what we noted is the underlying coupling between the myelination covariance and the functional synchronization for each paired region. We constructed these myelination covariant profiles based on the homophilic principles of nervous system organization, where the connections represent myelination covariant association that reflects the homogeneity of myelination growth between brain regions (the biological significance of covariant connectivity has been validated by previous studies19). The results, presented in the Section \"The dominant myelination-function coupling in the sensory cortex\" of the Results, confirm that myelination homogeneity supports the early formation and development of functional brain networks during infancy.\nTo further address your concern, we calculated the coupling between the functional connectivity profile and the direct spatial changes of myelination, as Myelin_FC. The Myelin_FC was calculated using the same method as MFC, adopting the linear regression model. Specifically, the predictor variable was the Myelin value of each vertex, and the dependent variable was the functional connectivity between a pair of vertices. The average value of individual Myelin_FC was significantly lower than that of MFC (t-value=101.86, p<0.0001, independent samples t-test; Figure R5a), and the Myelin_FC of all vertices across the whole brain was lower than the corresponding MFC. The increment of MFC compared with Myelin_FC could reach up to 438% (2nd-98th percentile range: 20%-270%; Figure R5b). The dominant regions with relative increase were distributed in the superior temporal gyrus, precentral sulcus, postcentral sulcus, cingulate sulcus, and higher visual cortex, whereas the primary sensorimotor cortex, parieto-occipital junction, and posterior cingulate cortex exhibited the lowest relative increase.\nFigure R5. Comparison between the Myelin_FC pattern and our MFC pattern. a Quantitative distribution of the two coupling indices. b Coupling increment map of MFC compared with Myelin_FC.\nThis suggested that myelination level itself contributed relatively little to the emergence of functional connectivity - information synchronization between spatially separated regions. In contrast, the myelination covariance between two regions acted as a marker of coordinated development of myelin and played a crucial role in enhancing the synchronization of functional activities between brain regions, which was consistent with the homophilic principles.\nWe wish the revised manuscript would be more convincing and readable than the previous version. If you believe our replies are not enough, please inform us to resolve them again.\nTo make it clear for readers, we have rewritten the Introduction (line 79 – line 132 in manuscript), and have expanded the Results by adding the Section “The dominant myelination-function coupling in the sensory cortex” (line 214 – line 217 in manuscript) and Supplementary Result 3 (see the Supplementary Information), which reads now as follows:\nIn the Introduction, “However, an important neural structural basis, myelination, in the infant brain is often overlooked. Myelination is known to be activity-dependent, thus coordinated developments with functional connectivity are expected 21–23. However, there is currently no previous study examining the direct, quantitative relationships between coordinated myelination and functional synchrony in the developing brain. Moreover, although white matter myelination has been more extensively studied, it remains unresolved whether it is the myelination of long-range tracts or the maturation of local gray matter that coordinates more closely with the emergent functional connectivity. These represent critical knowledge gaps in better understanding of infant structural-functional development.”\nIn the Results, “Furthermore, a direct coupling between the functional connectivity profile and the myelin map showed significantly lower coupling than our proposed MFC (see Supplementary Fig. 5 and Supplementary Result 3).”\nIn the Supplementary Result 3, “Another important question is whether the homogeneity of myelination between a pair of brain regions or vertices has greater functional significance than myelination itself in driving the emergence of functional connectivity. Following the method used to calculate MFC, we computed the coupling between myelination itself (as opposed to myelination covariance) and functional connectivity, with this relationship characterized by adjusted R² and designated as Myelin_FC. The average value of individual Myelin_FC was significantly lower than that of MFC (p<0.0001; Supplementary Fig. 5a), and the Myelin_FC of all vertices across the whole brain was lower than the corresponding MFC. The increment of MFC compared with Myelin_FC could reach up to 438% (2nd-98th percentile range: 20%-270%; Supplementary Fig. 5b). The dominant regions with relative increase were distributed in the superior temporal gyrus, precentral sulcus, postcentral sulcus, cingulate sulcus, and higher visual cortex, whereas the primary sensorimotor cortex, parieto-occipital junction, and posterior cingulate cortex exhibited the lowest relative increase.”\nComment 3. The interpretation of the T1w/T2w ratio as a measure of myelination is still debated. While more specific measures (e.g., R1 or MWF) are not available in the dHCP, I think this limitation needs to be explicitly acknowledged and discussed in more depth.\nResponse: We thank you for highlighting this important methodological consideration regarding the interpretation of the T1w/T2w ratio. We fully agree that while the T1w/T2w ratio is a widely used and sensitive proxy for cortical myelin content, it is not a specific measure and can be influenced by other tissue properties such as iron deposition and axon density.\nIn direct response to this comment, we have added a dedicated paragraph in the limitations of the Discussion to explicitly acknowledge and discuss this point in greater depth. The added text (line 570 - line 576 in manuscript):\n“Second, the T1w/T2w ratio was used to quantify myelin content in this study. Given its ease of calculation from commonly acquired images, T1w/T2w is usually used in large-scale datasets and/or early life samples where acquisition times are particularly restricted. It must be admitted that the validity of T1w/T2w as a marker of myelin remains debated72,73, although studies comparing T1w/T2w to R1 have validated its usability in adult cortex74 and infant white matter75.”\nComment 4. To my understanding, the authors assess covariance in functional activity and covariance in myelination. I find it difficult to follow the interpretation of this as “connectivity.” Could the authors clarify what is meant by “myelination connectivity”?\nResponse: Sorry for our unclear descriptions. \"Connectivity\" is a core noun in brain network research. Brain network analysis based on graph theory models the brain as a graph composed of nodes and edges: nodes generally represent brain regions or vertices, while edges denote connectivity in functional, anatomical, or other dimensions. For the functional connectivity (FC) network constructed in this study, nodes correspond to cortical vertices, and edges represent the synchronization of functional activities between vertex pairs, measured by Pearson correlation between the pairwise BOLD time series.\nRegarding the myelination covariance network, it refers to a brain network constructed by measuring correlations in structural or morphological features (e.g., cortical thickness, gray matter volume, here the myelination level is used) across different brain regions among a group of individuals. It reflects coordinated patterns of structural variation, potentially influenced by genetic factors, neurodevelopment, aging, or disease. Structural covariance network (SCN) provides insights into the large-scale organization of the brain’s structure, complementing functional and anatomical connectivity studies. Since the early 2000s, structural covariance network (SCN) analysis has developed rapidly as a complement to DTI tractography measures of structural connectivity20–23. Based on the known neurodevelopmental process of spatially organized, heterochronic maturation of the human cortex7,24, structural covariance analysis characterizes covariance by calculating the statistical correlations of structural indices between brain regions. The various methods for calculating covariant connectivity have been developed, primarily encompassing both population-level19 and individual-level dimensions, including the multi-metric-based MSN approach25, the multi-depth-based MPC method26, and the KL divergence-based MIND method27.\nIn this study, we integrated covariance at both the group-level and individual-level and proposed a dual-axis covariance framework (gMC and sMC) to investigate myelination-function coupling (MFC). Regarding interpreting the structural covariance as structural connectivity, the covariance essentially reflects the coordinated neurodevelopment across brain regions19. Therefore, myelination connectivity in the current study means the covariant level of myelination between two regions.\nTo help readers grasp the core idea of this research, we should provide a clearer explanation of myelination covariance. We now rewritten this section in the Introduction (line 109 - line 132 in manuscript):\n“Homophilic principles of nervous system organization state that brain regions with similar developmental trajectories will be similar also in terms of other properties34–36, hypothetically owing to shared genetic and/or convergent environmental effects. Based on the known neurodevelopmental process of spatially organized and heterochronic maturation of the human cortex35,37, structural covariant connectivity characterizes the coordination of neural development measured by statistical similarity in developmental trajectories or distributions of structural indicators between paired brain regions. Over the past two decades, various methods of structural covariance analysis have been developed 38–41, primarily encompassing both population-level and individual-level dimensions, and were used to study heterochronic network development42,43. Emerging evidence implicates gray matter myelination as a pivotal yet understudied mediator21,28,44,45. Cortical myelination is correlated with cortical development and intracortical circuit complexity46,47, and gray and white matter myelination are significantly correlated48. The first wave of myelination initiates prenatally, during which the primary sensory and motor cortices undergo myelination, and heteromodal association cortices do not become fully myelinated until the second decade, which is similar to the hierarchical development of cortical function 5,21,49. Additionally, the structural covariance connectome also has similar dense characters 50 with the functional connectome that even strong connectivity could be performed without direct fiber connectivity. Therefore, we held the hypothesis that the homophilic cortical growth measured by the myelination covariance during development could support the emergence of functional synchronization during infancy.”\nWe wish these replies have resolved your concern. If you believe our replies are not enough, please inform us to resolve them again.\nComment 5. The reported association between functional connectivity and myelination seems to be driven largely by neighboring vertices. Could their main observations simply be a consequence of the fact that both measures change smoothly in space?\nResponse: We thank you for raising this critical point regarding spatial autocorrelation, which is a key consideration in this surface-based coupling analysis. We agree that the smooth spatial variation of both myelination and functional connectivity maps could, in principle, drive spurious correlations between neighboring vertices. We have taken this concern seriously and performed multiple rigorous analyses to test and control for this effect.\nFirst, the proposed myelination function coupling (MFC) is the relationship between the myelination covariance and functional connectivity of the paired regions, not direct myelination, with our assumption that the homophilic development of brain structure will associate with the functional synchronization between two regions. Together with the MFC map in Figure 2b, we have also presented the significance of the MFC calculation model using spin permutation tests and most of the vertices showed significance (Supplementary Fig. 1b).\nSecond, as presented in the manuscript (Figure 4d and e), we explicitly regressed out the effect of spatial distance or removed the local connections from the myelination-function coupling (MFC) maps. These distance-controlled MFC maps retained a highly similar spatial pattern to the original maps (correlation r = 0.7 and r=0.88, Supplementary Fig. 7b), and the core findings regarding regional dominance (e.g., in primary sensory cortices) and system-level differences remained robust. Third, in the revised manuscript, we added more results with distance controlling, including the distance dependence analysis (see Supplementary Fig. 7c) and the age effects (see Supplementary Fig. 7d and e). These all demonstrate that the observed coupling is not a mere artifact of smooth spatial gradients.\nFurthermore, the distance dependence analysis provided the dependency of MFC strength as a function of cortical distance (with or without the distance control), by dividing the whole cortical vertices into several subsets according to the geodesic distance from local to remote (Figure 4a and b). If the coupling were purely driven by smooth spatial changes, we would expect a rapid, monotonous decay of MFC with distance. Instead, we observed a more complex, system-specific profile: a sharp initial drop followed by a stable, non-zero plateau for intermediate to long distances or even rise for specific regions.\nFinally, we wish to offer a conceptual reframing. While we have rigorously controlled the spatial confounding effect, we also posit that spatial proximity itself is a fundamental biological organizer in early brain development. Neighboring cortical regions are more likely to share developmental timing (heterochrony), genetic expression profiles, and exposure to similar microenvironmental cues.\nOur MFC metric captures how this coordinated microstructural maturation scaffolds the emergence of functional synchronization, starting from local neighborhoods and extending to distributed networks.\nComment 6. In the transcriptomic analyses, I was surprised that no genes related to myelination were identified. This absence seems important and deserves discussion.\nResponse: Thank you for this insightful comment. We agree that the absence of genes related to myelination in the transcriptomic analyses is a notable observation, and we appreciate the opportunity to discuss this finding in depth. Below is a detailed explanation.\nFirst, constrained by the limited number of postmortem tissue samples from the BrainSpan Atlas of the Developing Human Brain, we selected postmortem postnatal subject samples corresponding to 24 post-conceptual weeks (pcw) and 37 pcw. Samples with insufficient regional coverage were excluded under the assumption that the BrainSpan cohort and our study population share homogeneity. RNA-Seq expression averaged values were extracted from 11 brain regions, including the primary sensorimotor cortices (M1C, V1C, A1C, S1C) and other higher-order cortices (see Table R1 for details).\nFor each tissue sample, RNA-Seq expression levels were quantified as Reads Per Kilobase of transcript per Million mapped reads (RPKM). In the original manuscript, we calculated the Pearson correlation coefficients between RPKM and MFC across 11 brain regions. Then, 90 genes with significant correlations were identified using a threshold of p<0.001, all exhibiting a correlation coefficient greater than 0.84. In contrast, genes directly related to myelination (e.g., myelin basic protein (MBP), myelin-associated oligodendrocyte basic protein (MOBP), myelin-associated glycoprotein (MAG), and proteolipid protein 1 (PLP1)) showed either weaker correlations or no significant correlation with MFC (See the Figure R6) compared to these 90 genes, and were therefore not classified as significant.\nFigure R6. Pearson correlation between RPKM of genes directly related to myelination and MFC.\nTable R1. Brain regions of postmortem tissue samples from BrainSpan\nRegions Descriptions\nA1C primary auditory cortex\nDFC dorsolateral prefrontal cortex\nIPC posteroventral (inferior) parietal cortex\nITC inferolateral temporal cortex\nM1C primary motor cortex\nMFC anterior (rostral) cingulate (medial prefrontal) cortex\nOFC orbital frontal cortex\nS1C primary somatosensory cortex\nSTC posterior (caudal) superior temporal cortex\nV1C primary visual cortex\nVFC ventrolateral prefrontal cortex\nSecond, it is critical to clarify that the MFC is a connectome-level coupling trait, and the myelination covariance network captures coordinated myelination development between a pair of regions, rather than direct measures of myelin content variations in single region. The interregional myelination coordination is likely mediated by upstream regulatory pathways (e.g., synaptic pruning or neurotrophic signaling) rather than direct myelination-related genes. Therefore, it is reasonable that although no genes directly linked to myelination were included in the 90 significantly associated genes, our enrichment analysis results still indicated that these significant genes were enriched in gliogenesis - a biological process indirectly related to myelination. Our GO results (Original Figure 6a in the main text) highlight that the enriched pathways, including transepithelial transport and gliogenesis, correspond to processes that support communication between regions:\n(1) Transepithelial transport: This process mediates nutrient or metabolic waste exchange across epithelial barriers (e.g., blood-brain barrier (BBB)), which is critical for delivering metabolic resources to brain regions.\n(2) Gliogenesis: Gliogenesis is a core event in brain development. Among glial cells, astrocytes provide structural and functional support to neurons, oligodendrocytes serve as the source of myelination in the central nervous system (CNS), and microglia modulate learning and neurogenesis. By orchestrating the spatiotemporal dynamics of glial cells in the CNS, gliogenesis facilitates the coordinated structural and functional development of the brain during infancy.\nThird, drawing on previous studies28,29, the partial least squares (PLS) method is more suitable than direct Pearson correlation for handling the high-dimensionality and collinearity inherent in gene expression data. Therefore, in the revised manuscript, we improved the study methodology by supplementing PLS-based transcriptomic analysis, which has been substituted for the original analysis in the main text. Specifically, to explore the transcriptomic basis of myelination-function coupling development during neonates, we conducted PLS regression analysis to detect the predominant gene expression organization. The first component of PLS (PLS1) depicted a gene expression pattern that was mostly associated with MFC distributions, which explained 20.83% of the variance (p <0.0001; Figure R7a). Meanwhile, under the PLS1, significant spatial association between the MFC score and the gene expression score was observed (corr=0.91, p = 0.0001; Figure R7b). Furthermore, by utilizing a univariate one-sample Z test, we ranked these candidate genes based on normalized PLS1 weights, identifying 352 positively weighted genes (PLS1+: Z>3, p <0.001; Figure R7c) and 112 negatively weighted (PLS1-: Z<-3, p <0.001; Figure R7c).\nWe utilized the Metascape30 to align PLS1+ and PLS1- genes with the gene ontology (GO) terms. Notably, 66% of all the terms identified in our original results remained significantly enriched in the PLS-based transcriptomic analysis, for example, “transepithelial transport”, “gliogenesis”, “gland development”, and “tube morphogenesis”. The PLS1+ gene list exhibited significant enrichment in terms related to angiogenesis, transmembrane transport, neural regulation, and immune response (the top 20 GO enriched clusters shown in Figure R7d as orange bars, and Figure R7e), while the PLS1- gene were associated with neuronal morphogenesis and neural regulation (see Figure R7d, blue bars).\nIn more detail, for the GO enrichments of PLS1+ genes:\n(1) Side of membrane: Side of membrane refers to the polar region on the extracellular side of the cell membrane, where abundant signaling receptors are concentrated. High expression of related genes enhances the fidelity of signal transmission.\n(2) Tube morphogenesis: Angiogenesis and vascular development establish pathways for the transport of nutrients, oxygen, and metabolic waste, thereby maintaining the homeostasis of the brain microenvironment. The neural tube structure plays a vital developmental role in neural patterning31 and neuron differentiation32, High expression of related genes promotes the initial generation of distinct neuronal subtypes33，indicating the assembly of functional neuronal circuits。\n(3) Transepithelial transport: This process mediates nutrient or metabolic waste exchange across epithelial barriers (e.g., blood-brain barrier (BBB)), which is critical for delivering metabolic resources to brain regions.\n(4) Gliogenesis: Gliogenesis is the biological process by which neural progenitor cells proliferate and differentiate into glial cells, which is crucial during brain development for constructing the brain structure34. Glial cells include three main types in the central nervous system (CNS): astrocytes, oligodendrocytes, and microglia35. Astrocytes are involved in synaptic function, oligodendrocytes participate in myelin sheath formation, and microglia modulate learning and neurogenesis36–38. Generally, gliogenesis follows neurogenesis with relatively little temporal overlap39. By orchestrating the spatiotemporal dynamics of glial cells in the CNS, gliogenesis facilitates the coordinated structural and functional development of the brain during infancy.\nFor the GO enrichments of PLS1- genes, genes annotated to neural structural components (e.g., “synaptic membrane”, “neuron projection cytoplasm”, and “axon”) were negatively correlated with MFC. This finding indicates that increased MFC corresponds to low expression of these genes, which supports synaptic pruning and neural structural maturation, and consequently contributes to the progressive stabilization of functional connectivity pathways.\nFigure R7. Transcriptomic decoding of MFC. a Variance explained for the first 10 components derived from the PLS regression analysis. b Pearson’s correlation between MFC score and gene expression score (r=0.91, p=0.0001). 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.\nWe wish these replies have resolved your concern. If you believe our replies are not enough, please inform us to resolve them again.\nTo make it clear for readers, we added the following content to the Results (section “The transcriptomic association of myelination-function coupling”, line 344 – line 364 in manuscript), Discussion (line 500 – line 532 in manuscript), and Methods (section “Spatial association with transcriptional gene expression”, line 765 – line 788 in manuscript):\nIn the Results, “To explore the transcriptomic basis of the myelination-function coupling during neonates, this study utilized the BrainSpan Atlas of the Developing Human Brain, aged from 24 post-conceptual weeks (pcw) to 37 pcw with 11 regions. The first component of partial least squares (PLS) depicted gene expression pattern that was mostly associated with MFC distributions, which explained 20.83% of the variance (p <0.0001; Fig. 6a). Meanwhile, under the PLS1 (first component of PLS), significant spatial association between the MFC score and the gene expression score was observed (corr=0.91, p = 0.0001; Fig. 6b). Furthermore, by utilizing a univariate one-sample Z test, we ranked these candidate genes based on normalized PLS1 weights, identifying 352 positively weighted genes (PLS1+: Z>3, p <0.001; Fig. 6c) and 112 negatively (PLS1-: Z<-3, p <0.001; Fig. 6c). We utilized the Metascape55 to align PLS1+ and PLS1- genes with the gene ontology (GO) terms, including biological process (BP), molecular function (MF) and cellular component (CC). The PLS1+ gene list exhibited significant enrichment in terms related to angiogenesis, transmembrane transport, neural regulation and immune response (the top 20 GO enriched clusters shown in Fig. 6d, orange bars, and Fig. 6e), such as “side of membrane”, “tube morphogenesis”, “transepithelial transport” and “gliogenesis”. The PLS1- genes were associated with neuronal morphogenesis and neural regulation (see the Fig. 6d, blue bars), such as “synaptic membrane”, “neuron projection cytoplasm” and “axon”.”\nIn the Discussion, “In addition, the transcriptomic associations identified in this study delineate a potential molecular blueprint underlying the observed myelination-function coupling, offering mechanistic insights into how coordinated myelination may drive the initial assembly of functional pathways. Our gene co-expression analysis revealed that regions with high MFC (PLS1+) are characterized by the upregulation of genes significantly enriched in biological processes crucial for constructing and maintaining the neurovascular coupling and the glial environment. Specifically, enrichments in “transepithelial transport” and “side of membrane” functions suggest the maturation of the blood-brain barrier (BBB) 65,66 and the cell membrane’s signaling apparatus, which are vital for maintaining ionic homeostasis and facilitating efficient neurotransmission. Crucially, the significant enrichment of “gliogenesis” terms, particularly those related to astrocyte and oligodendrocyte precursor development, bridges these vascular and neuronal components. Astrocytes are central players in neurovascular coupling 9,15, dynamically regulating local blood flow in response to neuronal activity, while oligodendrocytes are responsible for myelination itself and microglia modulate learning and neurogenesis67–69. In contrast, the genes whose expression was negatively associated with MFC (PLS1−) were enriched for components intrinsic to the post-mitotic neuronal structure and its refinement, such as “synaptic membrane”, “neuron projection cytoplasm” and “axon”. This pattern may reflect a complementary, and perhaps temporally later, phase of circuit maturation. The relative downregulation of genes involved in axonal growth and synaptic adhesion in high MFC regions could indicate a state of increased circuit stability or a shift towards activity-dependent pruning.23 In this view, once the foundational neurovascular-glia unit is established via PLS1+ gene activity, securing basic metabolic and signaling efficiency, the focus of development may shift towards the activity-dependent refinement of neuronal connectivity. The lower expression of PLS1- genes might therefore signify a transition from exuberant growth to the selective stabilization of the most functionally relevant synapses and projections within the co-varying regions. This molecular model provides a testable hypothesis for how myelination covariance scaffolds the sequential stages of functional pathway formation from initial synchronization to refined efficiency.”\nIn the Methods, “To elucidate the underlying biological mechanisms of myelination-function coupling during neonates, this study utilized transcriptomic data from donors aged 24 post-conceptual weeks (pcw) and 37 pcw obtained from the BrainSpan Atlas (https://www.brainspan.org/static/home) of the Developing Human Brain. Samples with insufficient regional coverage were excluded under the assumption that the BrainSpan cohort and our study population share homogeneity. RNA-Seq expression averaged values were extracted from 11 brain regions, including the primary sensorimotor cortices (M1C, V1C, A1C, S1C) and other higher-order cortices (see the Supplementary Table 7 for details). Then, PLS (partial least squares) regression was employed to probe the potential relationship between MFC and transcriptional levels for all genes. PLS1(first component of PLS) depicted the most correlated cortical expression map with respect to the MFC. Additionally, we performed a nonparametric permutation test by randomly shuffling the orders (participant identities) of regional MFC 5,000 times to assess the statistical significance of PLS1. The bootstrapping technique was performed with 5000 bootstrap replicates to evaluate the contribution of each gene to the PLS1. The Z-score per gene was determined by dividing its weight by the corresponding bootstrap standard error, based on which all genes were ranked. Then, genes exhibiting Z>3 and p <0.001 were assigned to the PLS1+ gene list, while those with Z<-3 and p <0.001 were allocated to the PLS1- gene list. Subsequently, the PLS1+ and PLS1− gene lists were separately subjected to the Metascape (https://metascape.org/gp-/index.html) to identify enriched GO terms, including biological process (BP), molecular function (MF), and cellular component (CC).”\nComment 7. For the behavioral associations, it would be interesting to examine whether networks predict their expected functional domains (e.g., whether motor system MFC predicts motor ability).\nResponse: We are grateful to you for this excellent and insightful suggestion. We agree that examining whether the MFC within specific functional networks preferentially predicts corresponding behavioral domains is a logical and powerful extension of this work, as it directly tests the functional specificity of our measure.\nFirst of all, establishing robust associations between the neonatal brain imaging features and the behavior scores at about 18 months of age is inherently challenging, since the neonatal brain undergoes fast growth and complicated interactions with the environment. As largely reported 74–77, these brain-behavior relationships are usually around 0.1~0.3. Notably, our PLS analysis, as presented in Figure R8a (enhanced from the Fig. 6 concerning behavior analysis in our original manuscript), captured the significant relationship between the proposed MFC of the whole brain and the behavior (p <0.01).\nSecond, and most critically, we directly followed your suggestion to perform a system-specific analysis. While a single model using the mean MFC of each system did not yield significant associations, a more refined, connection-distance-dependent way revealed a precise and theoretically meaningful association. Specifically, we found that only remote distance (long-range connection) MFC within the sensorimotor network (SMN) showed a significant association with 18-month motor capability (adjusted R2=0.045, p <0.05; Figure R8c). This was the sole significant finding across all functional systems and behavioral domains.\nFigure R8. Relationship between MFC and the behavior scores at 18 months of age. a Significant association between behavioral and MFC composite scores under the PLS1 (adjusted R2=0.026; p =0.006), with the behavioral and MFC loadings. b Significant association between behavioral and DFC composite scores in the first latent component (adjusted R2=0.047; p =0.001), with the behavioral and DFC loadings. c Significant association between the motor domain and MFC of motor network composite score in the first latent component (adjusted R2=0.045; p =0.041), with the domain loading and the MFC loadings. d The monotonic relationship between the R2 and the distance (r=0.773, p=0.015), together with the MFC loadings for each distance.\nThird, we interpret this result as strong evidence for domain-specific relevance. It indicated that the coupling between myelination covariance and function is not uniform. The integrity of long-range integrative connections within the neonatal motor system specifically scaffolds the subsequent development of motor skills. This finding significantly strengthens the biological relevance of our MFC framework by moving from a generalized brain-behavior correlation to a domain and spatial-specific prediction.\nTo make it clear for readers, we have added details to the Results (section “The behavioral association of myelination-function coupling”, line 366 – line 383 in manuscript), Discussion (line 533 – line 559 in manuscript):\nIn the Results, “Behavioral assessments from the Bayley III scale for dHCP subjects at approximately 18 months of age were associated with the individual MFC values in this study. Non-parametric permutation analysis (5000 times permutation) for null-hypothesis testing revealed a statistically significant population correspondence between the regional MFC and behavioral profiles. Whole-brain MFC and behavioral scores for PLS1 showed a significant correlation (adjusted R2=0.026, p =0.006; Fig. 7a), with the negative behavior loadings and positive image loadings. Similarly, under the PLS1, whole-brain functional connectivity strength (FCS) and behavioral scores also showed significance (adjusted R2=0.047, p =0.001; Fig. 7b), with corresponding behavior and imaging loadings. For functional system-level PLS analysis, remote-distance sensorimotor MFC was significantly associated with the motor behavioral domain (adjusted R2=0.045, p =0.041; Fig. 7c). Notably, for the 9 distance subsets (see Fig. 4a), the distance-specific MFC-behavior relationships (adjusted R2 from the first latent variable of PLS) exhibited increasing association with the connection distance, with adjusted R2 values increasing from 0.020 for the local connections to 0.092 for the remote connections (Fig. 7d and Supplementary Fig. 10).”\nIn the Discussion, “Another notable validation of our MFC framework lies in its significant association with later behavioral outcomes at about 18 months of age. These associations are not merely global but exhibit a striking functional and anatomical specificity that aligns with the proposed developmental mechanism. At the whole-brain level, MFC showed a comparatively lower association with behavioral outcomes than FCS. This pattern is informative rather than contradictory, as it likely reflects the distinct neurobiological processes captured by each metric and the developmental stage of our cohort. FCS predominantly measures the magnitude of spontaneous neural co-activation, and MFC quantifies the alignment between myelination covariance and this functional co-activation. In the neonatal period, this alignment may be a more nascent and regionally specific process, and only index explains a significant yet partial portion of the variance in functional connectivity. Functional networks are supported by a multiplex of structural substrates that operate across different spatial and temporal scales. Our study highlights intracortical myelination covariance as a critical, early-developing component of this multiplex, particularly relevant during the neonatal period of profound white matter immaturity. Furthermore, the strength of the MFC-behavior association exhibited a clear distance-dependent gradient, increasing monotonically from the shortest to the longest connections. This observation indicates that long-distance integration, which involves the synchronization of distant brain regions, plays a crucial role in generating behaviorally relevant connections.9,70 Long-distance coupling supports the notion that the integration of more distant cortical areas contributes to the development of complex cognitive functions and behaviors. These findings emphasize that long-range connections, despite their relatively later myelination, are integral to shaping the brain's ability to coordinate complex behaviors 8,56,70. ”\nComment 8. Please include p-values in the main text.\nResponse: Thanks for your careful checks. We have added p-values to the main text. And here we did not list the changes, but marked them in red in the revised manuscript.\nComment 9. Could the authors explain why T1w/T2w measures were z-scored? I am concerned that this procedure might obscure meaningful inter-subject differences, such as those related to age.\nResponse: Sorry for not providing a clear explanation in the original manuscript. Actually, in the calculation of covariance, z-normalization was performed on T1w/T2w measures within each individual to eliminate baseline differences in individual data, thereby making the relative patterns of brain regions comparable across different individuals.\nThe advancements driven by researchers worldwide have led to the development of multiple methodologies for individual structural covariance, including the multi-metric-based MSN approach25, the multi-depth-based MPC method26, and the KL divergence-based MIND method27. As presented in our original manuscript, MFC derived from the KL divergence-based approach at the atlas level was highly similar with the MFC computed via our proposed method (corr=0.81, p<0.0001; Supplementary Fig. 2).\nFor the calculation of sMC, T1w/T2w measures were first z-scored within each individual and then subjected to dot product. This design offers two critical benefits:\n(1) The sign of the dot product can indicate the presence or absence of developmental coordination. After z-scoring, a positive dot product indicates that the myelination levels of the paired vertices are coordinated—both being higher or lower than the average level—whereas a negative dot product signifies uncoordinated myelination, with one vertex above the average level and the other below.\n(2) After z-scoring, the dot product values are comparable across individuals, eliminating the confounding effect of individual baseline variability on the relative differences between paired vertices. Higher values correspond to stronger intra-subject coordination of the paired vertices, while lower values reflect weaker coordination.\nTherefore, z-normalization is necessary. Moreover, following this appropriate transformation, the age effect captures age-related changes in the intra-subject coordination pattern of each vertex relative to other vertices—an observation that is more meaningful than age-related variations in the product of absolute values.\nFor the calculation of group-level myelination covariance (gMC), the z-score operation was performed to align with the sMC calculation procedure. This operation did not induce substantial changes in the MFC results. At both the vertex-level and global-level, a high degree of similarity was observed when comparing MFC derived from gMC with and without the z-normalization operation (vertex-level: corr=0.95, p<0.0001, Figure R9a; global-level: corr=0.98, p<0.0001, Figure R9b).\nFurthermore, when modeling the age effect of MFC using generalized additive model (GAM), the age effect of the MFC derived from gMC without z-normalization was highly consistent with the original findings (Figure R9c and R9d). Dice coefficient between the significantly developing brain region map and the original map was as high as 0.91 (Figure R9c and Figure 3b).\nFigure R9. Comparison between the original MFC and the MFC computed from gMC without z-score. a. Pearson correlation between the original group averaged MFC map and the group averaged MFC without z-score at the vertex level (r=0.95, p<0.0001). b. Pearson correlation between the original individual MFC and the individual MFC without z-score at the global level (r=0.98, p<0.0001). c. Growth-rate coefficient of MFC without z-score derived from vertex-wise GAM analysis with p <0.05. d. Gobal-level growth rate of MFC without z-score.\nWe wish these replies could resolve your concerns. If you believe this is not enough, please inform us, and we will try our best to resolve this issue again.\nReviewer #2 (Remarks on code availability):\nI did not try to run the code, I just looked at it briefly. It looks well organized and easy to follow. However, it seems not all analyses are included at the moment? E.g. I could not find code to run the behavioral associations.\nResponse: Thanks. Basically, we would like to build and share the entire code including the whole analysis of our work on GitHub. Now, in the updated codes folder, we uploaded the toolbox as well as our scripts, together with their reference citation and link address in the CODE AVAILABILITY. For the behavior analysis with PLS regression, we created the analysis based on a public toolbox (https://github.com/danizoeller/myPLS ).\nReviewer #3 (Remarks to the Author):\nThis paper investigates the relationship between myelination and functional connectivity (FC) during infancy using fMRI data. While the topic is of interest, several significant methodological and analytical concerns need to be addressed.\nResponse: Thanks for your insightful comments. We have carefully considered these comments and accordingly revised our work. In this reply letter, we have also addressed all these comments one by one. Please kindly refer to these replies for modifications made in our revised manuscript.\nComment 1. Head Movement Artifacts in Infant fMRI: Infants are prone to significant movement, and the manuscript does not provide a robust description of the head movement correction pipeline. Large head movements can introduce spurious functional correlations that could confound the central findings.\nResponse: Thanks for your kind comment about the head motion. It is indeed well documented that fMRI analysis is very sensitive to subject head motion79. Since the data is publicly released data, we did not fully preprocess the data from the raw data, but started from the minimally preprocessed data. That’s the reason why we did not write sentences about the head motion in the previous version. Firstly, the minimal preprocessing pipeline has considered to carefully control the potential head motion effects in fMRI images, including the motion and distortion correction (MCDC) with FSL EDDY, the denoising based on spatial independent component analysis (sICA) and nuisance regression with motion parameters, temporally band-pass filtered (0.01-0.09Hz) and spatial registration80. The MCDC with FSL EDDY is a novel way to control the intra-volume movement artefacts81. Spatial independent component analysis (sICA) is also performed to remove the residual motion related artefacts, as well as the other additional structured noise artefacts which are unrelated to head motion. The FIX labeled artefact ICs and motion parameters were simultaneously regressed from the motion and distortion corrected functional time-series.\nSecondly, in order to estimate the residual relationship between subject head motion and functional connectivity after the preprocessing, we computed QC-FC correlations (mFC - FC), and the distribution of correlations was zero-centered with no significance. The age (PMA) effect of the residual mFDs was also not significant (p=0.435). Thirdly, in the following group analysis, the head motion parameter (the residual mFD) was considered as one of the covariates.\nTo make it clear for readers, we have also added the following content to the Methods (section “MRI Data Preprocessing”, line 624 – line 633 in manuscript):\n“For the rs-fMRI data, the minimally pre-processed 4D functional volume images which were motion and distortion corrected (MCDC) with FSL EDDY, denoised based on spatial independent component analysis (sICA) and nuisance regression, temporally band-pass filtered (0.01-0.09Hz) and spatial registration. Even though the potential head motion effects on fMRI signals had been carefully controlled in MCDC, sICA and motion parameters regression, we further computed the relationship between mean relative RMS displacement and functional connectivity (QC-FC correlations)82 and the distribution of QC-FC correlations was zero-centered with no significance. The residual mFDs were further examined for age effects and nothing significant was found. For the registration challenge since…”\nComment 2. Exclusion of Subcortical Structures: Early brain development is characterized by the maturation of deep-seated subcortical structures and their connections to the cortex. For example, the early development of connections between the thalamus and cortex might be a covariant cause for both functional connectivity and myelination development. The observed covariance between myelination and FC could be driven by these subcortical connections, and their exclusion from the analysis could lead to an incomplete and potentially misleading conclusion.\nResponse: We are grateful to you for this constructive suggestion, and we definitely agree that the paper would benefit from the consideration of the influence of subcortical structures. Subcortical nuclei develop early in the developmental trajectory83,84 and form extensive projection connections with cortical regions85,86. We supplemented a dedicated analysis of subcortical nuclei to investigate their impact on the coupling of myelination covariance with functional connectivity.\nFirst, the MFC calculation results with and without the subcortical nuclei regions were compared. We analyzed the results of two atlases: the UNC infant AAL atlas (comprising 8 subcortical structures and 82 cortical structures) and the INFANT atlas15 (comprising 12 subcortical structures and 211 cortical structures). Myelination map and BOLD time series from either all brain structures, including subcortical regions, or only cortical regions were utilized to construct the two axes of MCN (Myelination Covariant Network, gMC and sMC) and FCN (Functional Connectivity Network). While the MFC incorporating subcortical structures shared a highly similar spatial distribution pattern with that derived from only cortical region data (AAL atlas: corr=0.98, p<0.0001, Figure R10a; INFANT atlas: corr=0.98, p<0.0001, Figure R11a), it is disproportionately diminished in every cortical region. Specifically, the most prominent decrement occurred in the dorsolateral prefrontal cortex and posterior cingulate cortex (Figure R10a, cortical coupling decrement map, 2nd-98th percentiles range: 5%~28%). This indicated that the coupling of myelination covariant connectivity between cortex and nuclei with functional connectivity between cortex and nuclei was lower than that at the cortical level, which was also confirmed at the global level (Figure R10c and Figure R10c): the relationship of two single-axis myelination covariant matrices with functional connectivity between the cortex and subcortical nuclei were lower than that observed within the cortex (AAL atlas: R2=0.21, beta=0.02, p<0.0001 for sMC,: R2=0.007, beta=-0.05, p<0.05 for gMC, Figure R10c; INFANT atlas: R2=0.089, beta=0.01, p<0.0001 for sMC, R2=0.11, beta=0.16, p<0.0001 for gMC, Figure R11c).\nSecond, we conducted another analysis about the impact of subcortical regions on our proposed MFC. The BOLD time series of all subcortical structures was included as a covariate to remove their neural signals from the BOLD time series of cortical regions, using partial correlations of functional activity signals between cortical regions. After this signal removal, the MFC calculation was performed. Notably, the updated coupling of the cortical functional connectivity (with subcortical structural influences excluded) and myelination covariant connectivity remained essentially unchanged compared to the original MFC results (AAL atlas: corr=0.99, p<0.0001, Figure R10d; INFANT atlas: corr=0.99, p<0.0001, Figure R11b).\nFigure R10. Comparison of MFC with and without the inclusion of subcortical nuclei at the AAL atlas. a Averaged map of individual MFC with and without inclusion of subcortical nuclei (SUB), and cortical decrement map of MFC-with-SUB compared with MFC-without-SUB. b Spatial association between MFC-with-SUB and MFC-without-SUB (r=0.98, p<0.0001). c Linear regression between mFC and gMC or sMC for the connections between the cortex and SUB. d Spatial association between MFC-without-SUB and MFC derived from FC with regressing out averaged BOLD time series of SUB (r=0.99, p<0.0001).\nFigure R11. The comparison of MFC with and without the inclusion of subcortical nuclei at INFANT atlas. a Pearson correlation between MFC-with-SUB and MFC-without-SUB (r=0.98, p<0.0001); b Pearson correlation between MFC-without-SUB and MFC derived from FC with regressing averaged BOLD time series of SUB (r=0.99, p<0.0001). c Linear regression between mFC and gMC or sMC for the connections between the cortex and SUB.\nOur findings indicate that subcortical structures have no significant impact on the coupling of cortical myelination covariance with functional connectivity during the early postnatal period. This observation could be well-supported by converging evidence from multiple domains:\nFirst, beyond differing timelines, the nature of maturation itself differs. During the late prenatal and early postnatal period, the cerebral cortex undergoes explosive, experience-expectant synaptogenesis and dendritic arborization, processes tightly coupled to the initial wave of intracortical myelination studied here. In contrast, key subcortical structures follow a more prenatal-centric schedule for core circuit formation and myelination, with much of their major structural scaffolding established before term-equivalent age.87 This fundamental developmental uncoupling means that the dynamic, regressive, and experience-dependent microstructural changes dominating the neonatal cortex are not occurring in synchrony with the relatively more stabilized subcortical maturation.\nSecond, empirical studies of neonatal and infant functional networks consistently report that the strongest and most rapidly consolidating functional connections are within the cortex, particularly in primary sensory and motor systems. While functional connections between the cortex and thalamus exist, they exhibit different developmental trajectories and lower strength compared to intra-cortical connections at this age. This supports the notion that the emergent functional architecture captured by resting-state fMRI in early infancy is primarily driven by and reflective of cortico-cortical coordination.\nThird, from a signal and modeling perspective, our dual-axis framework quantifies covariance in cortical myelination profiles. The spatially diffuse and polysynaptic nature of subcortical influence on widespread cortical territories may not be effectively captured by a pairwise covariance model designed to assess homophilic, large-scale cortical gradients. The observed lack of significant impact may reflect that subcortical modulation operates through a different mechanism that is not well-modeled as a direct contributor to the covariance of cortical microstates.\nTo make it clear for the reviewer and readers, we have added the following paragraph in the Results (section “The dominant myelination function coupling in the sensory cortex”, line 211 – line 214 in manuscript) and Supplementary Result 2 (see the Supplementary Information):\nIn the Results, “Considering the subcortical regions, the MFCs were highly consistent (corr=0.98, p<0.001) between that eighter with or without the subcortical regions in the calculating framework (Supplementary Fig. 4 and Supplementary Result 2).”\nIn the Supplementary Result 2, “To enhance the completeness of our research, we also explored the effects of early-developing subcortical structures on the MFC. First, the MFC calculation results with and without the subcortical nuclei (SUB) regions were compared. For the UNC infant AAL atlas, myelination map and BOLD time series of all 90 brain structures (including 8 subcortical structures and 82 cortical structures) were used to construct the two-axes MCN (Myelination Covariant Network, gMC and sMC) and FCN (Functional Connectivity Network). Compared with the MFC derived from data of only 82 cortical regions, the MFC incorporating subcortical structures exhibited a highly similar spatial distribution pattern (corr=0.98, p<0.0001, Supplementary Fig. 4a) but was reduced disproportionately in all cortical regions, with the most prominent decrement occurring in the dorsolateral prefrontal cortex and posterior cingulate cortex (Supplementary Fig. 4b, cortical coupling decrement map, 2nd-98th percentiles range: 5%~28%). This indicated that the coupling of myelin covariance connectivity between cortex and nuclei with functional connectivity between cortex and nuclei was lower than that at the cortical level, which was also confirmed at the global level (Supplementary Fig. 4c): the relationship of two single-axis myelination covariant matrices with functional connectivity between the cortex and subcortical nuclei were lower than that observed within the cortex (R2=0.21, beta=0.02, p<0.0001 for sMC; R2=0.007, beta=-0.05, p<0.05 for gMC). Furthermore, the mean BOLD time series of all subcortical structures was included as a covariate to remove their neural signals from the BOLD time series of cortical regions, using partial correlations of functional activity signals between cortical regions. Notably, the updated coupling of the cortical functional connectivity (with subcortical structural influences excluded) and myelination covariant connectivity remained essentially unchanged compared to the original MFC results (corr=0.99, p<0.0001, Supplementary Fig. 4d).”\nComment 3. Novelty and Differentiation is limited: The reported patterns of \"coupled regions\" in the results appear to be broadly consistent with previous literature that sensory systems are developed earlier. It is difficult to conclude that the findings are uniquely due to myelination or other well-known developmental factors. The analysis would be significantly strengthened by leveraging the Diffusion Tensor Imaging (DTI) data available in the Developing Human Connectome Project (dHCP) dataset to also model white matter tract development and its relationship with the observed functional changes.\nResponse: We sincerely thank you for this incisive critique, which compels us to articulate the conceptual novelty of our work with greater clarity and to strengthen it through the suggested integration of white matter data. We address your concern in turn below.\nFirst, we agree that the earlier maturation of sensory systems is a well-established principle. The novelty of our study lies not in rediscovering this hierarchy, but in proposing and validating a novel, mechanistic framework to explain how functional networks can emerge despite the profound immaturity of the white matter infrastructure that is classically assumed to support them. While previous studies have investigated the coupling of gray matter structural covariance and cortical function using different structural covariant connections (e.g., SCN19, MPC26, MSN25, and MIND27), no previous work, to our knowledge, has quantified the association between gray matter myelination covariance and cortical function in infancy, nor has a study integrated covariant connections at the population level and individual level.\nOur work is grounded in two key observations that create a compelling neurodevelopmental paradox: (1) myelination levels of neonatal white matter tracts are far-from-mature (30% of the adult level), yet (2) neonates exhibit resting-state functional networks sharing 80% spatial similarity with adult counterparts. This suggests that white matter-based structure-function coupling is likely insufficient to fully explain the underlying dynamics of the developing infant brain. We therefore introduced the dual-axis myelination covariance framework, hypothesizing that the synchronized microstructural maturation of cortical regions - a reflection of shared developmental trajectories under the “homophilic principle”- provides a complementary, intrinsic blueprint for initial functional synchronization. To our knowledge, this is the first study to quantify the association between gray matter myelination covariance and functional connectivity in infancy, and the first to integrate both group-level and subject-level covariance within a unified model.\nSecond, our findings extend beyond confirming a known sensory-to-association gradient. They provide new, quantitative insights into the spatial, temporal, and behavioral dimensions of this early structure-function relationship:\n(1) Spatial specificity. We map the precise cortical landscape of myelination-function coupling (MFC), identifying primary sensory/motor and insular regions as hubs.\n(2) Developmental dynamics. We show that MFC undergoes rapid, regionally heterogeneous growth modulated by postnatal age, the birth transition, and connection distance, capturing dynamic maturation in the first weeks of life.\n(3) Multimodal validation. We demonstrate that neonatal MFC patterns are associated with relevant gene expression profiles and, crucially, show that long-range MFC growth predicts better 18-month behavioral outcomes.\nThis links our microstructural metric to molecular mechanisms and long-term behavioral development, moving from correlation to mechanistic and functional relevance.\nThird, we fully agree with your suggestion that integrating white matter tract data is essential to differentiate our MFC from a mere reflection of underlying tract development. In direct response, we have conducted comprehensive new analyses using dHCP DTI data (consistent with our responses to Comments #1 of Reviewer #2).\nWe constructed a white matter structural connectome, derived white matter-based structure function coupling (SFC) index, and performed comparative and integrative analyses (see Figure R3). The key results are:\n(1) Spatial and quantitative distinction: The SFC map was more limited in extent and significantly weaker in global strength than the MFC map (mean SFC = 0.09 vs. MFC = 0.22).\n(2) Unique predictive or informative value: For most of the cortex, adding white matter connectivity to our model did not substantially improve the prediction of functional connectivity beyond the variance explained by gray matter myelination covariance alone.\nThese results, to be included in the revised manuscript, provide strong empirical evidence that MFC captures distinct variance in neonatal functional organization that is not redundant with, and is more spatially extensive than, contemporaneous white matter tract strength. This directly supports our central topic that coordinated gray myelination maturation provides a primary and distinct scaffold for early functional network emergence. We have added details in the Results (section “The dominant myelination function coupling in the sensory cortex”, line 206 – line 211 in manuscript) and Supplementary Result 1 (see the Supplementary Information):\nIn the Results, “Regarding the white matter effect, we calculated the white matter-based structure-function coupling (SFC) map and it showed a similar spatial pattern (corr=0.48, p<0.001) with our proposed MFC map, but about eight times lower level than MFC (see Supplementary Fig. 3 and Supplementary Results 1). Including structure connectivity (SC) profiles in our MFC framework only gains limited benefits to the coupling level.”\nIn the Supplementary Result 1, “Having established that structural covariance of gray matter myelination is associated with the emergence of functional connectivity, it remained possible that white matter fiber connections also play a role therein. Therefore, we incorporated DTI data for fiber tract reconstruction, calculated the fiber tract connectivity density between brain regions to construct the structural connectivity network (SCN) based on the UNC infant AAL atlas (see Supplementary Method 1 for details). Subsequently, the fully weighted SCN was used to model the coupling relationship between white matter tracts and functional connectivity. The coupling was calculated using the same method as MFC, adopting a linear regression model. Averaged map of atlas-level individual MFC and SFC were illustrated (Supplementary Figure 3a and b), showing that a moderate level of correlation between them (corr=0.48, p<0.0001; Supplementary Fig. 3c). MFC map showed much higher level than SFC across the whole brain (the mean MFC is 0.22, and the mean SFC is 0.09), except for the cuneate cortex, with the peak increment up to 9 times (Supplementary Fig. 3d, the coupling increment map of MFC over SFC). The dominant regions of SFC were only distributed in the medial cortex, including the primary visual cortex and primary motor cortex. Furthermore, we introduced structural connectivity (SC) as an additional predictive variable into the MFC dual-axis framework. For the majority of brain regions, the introduction of SC profile barely improved the prediction strength for functional connectivity, particularly in the dorsolateral prefrontal cortex (Supplementary Fig. 3e, MFC map with SC included). Regions with an increment exceeding 20% were concentrated in the posterior cingulate cortex, calcarine cortex, and cuneus (Supplementary Fig. 3f, Coupling increment map of SC-included MFC over MFC without SC variable).”\nAdditionally, we have rewritten the Introduction (line 79 – line 132 in manuscript) to more clearly present the research hypothesis, and to better communicate the conceptual novelty of this integrated neonatal myelination-function coupling study as,\n“However, an important neural structural basis, myelination, in the infant brain is often overlooked. Myelination is known to be activity-dependent, thus coordinated developments with functional connectivity are expected 38,73,87. However, there is currently no previous study examining the direct, quantitative relationships between coordinated myelination and functional synchrony in the developing brain. Moreover, although white matter myelination has been more extensively studied, it remains unresolved whether it is the myelination of long-range tracts or the maturation of local gray matter that coordinates more closely with the emergent functional connectivity. These represent critical knowledge gaps in better understanding of infant structural-functional development.\nAs extensively reported, myelination is immature in the newborn stage. However, despite far-from-mature myelination levels in white matter tracts 24,25 (30% of the adult level), neonates exhibit resting-state functional networks sharing 80% spatial similarity with adult counterparts26,27. This suggests that white matter-based structure-function coupling is likely insufficient to fully explain the underlying dynamics of the developing infant brain. Therefore, we hypothesize that the formation of functional connectivity does not solely rely on white matter fiber connections. Instead, we propose that coordinated growth of local gray matter myelination is more closely related to the emergence of initial functional connectivity. The brain wiring, including synaptic formation and pruning within the signaling-integrated white matter, is more potently modulated by the distributed gray matter systems 28–32. Notably, gray matter, which primarily consists of neuronal cell bodies, develops earlier in the fetus than white matter, and undergoes rapid growth in late gestation and continues to develop after birth.5–7 Unlike white matter myelination that peaks postnatally, gray matter oligodendrogenesis commences as early as 28 weeks of gestational age, synchronizing with the emergence of local field potential oscillations1,28,29,33,34. Such findings imply that early functional architectures may arise from the cumulative cortical growth.\nHomophilic principles of nervous system organization state that brain regions with similar developmental trajectories will be similar also in terms of other properties34–36, hypothetically owing to shared genetic and/or convergent environmental effects. Based on the known neurodevelopmental process of spatially organized and heterochronic maturation of the human cortex35,37, structural covariant connectivity characterizes the coordination of neural development measured by statistical similarity in developmental trajectories or distributions of structural indicators between paired brain regions. Over the past two decades, various methods of structural covariance analysis have been developed 38–41, primarily encompassing both population-level and individual-level dimensions, and were used to study heterochronic network development42,43. Emerging evidence implicates gray matter myelination as a pivotal yet understudied mediator21,28,44,45. Cortical myelination is correlated with cortical development and intracortical circuit complexity46,47, and gray and white matter myelination are significantly correlated48. The first wave of myelination initiates prenatally, during which the primary sensory and motor cortices undergo myelination, and heteromodal association cortices do not become fully myelinated until the second decade, which is similar to the hierarchical development of cortical function 5,21,49. Additionally, the structural covariance connectome also has similar dense characters 50 with the functional connectome that even strong connectivity could be performed without direct fiber connectivity. Therefore, we held the hypothesis that the homophilic cortical growth measured by the myelination covariance during development could support the emergence of functional synchronization during infancy.”\nWe hope this revision addresses your concerns and clarifies the scientific question and motivation of our current work. If more modifications are still needed, please contact us to resolve this again.\nComment 4. The relationship between MFC and Distance: The manuscript notes a strong presence of MFC in short-distance connections. The paper does not adequately explore whether this is due to a lack of myelination in white matter, which would imply that myelination is more critical for long-range FC development. It is highly suggested that distance be defined based on fiber length or myelination level, rather than a simple geodesic distance, to provide more meaningful insights into this relationship.\nResponse: We sincerely thank you for this insightful suggestion. We completely agree that defining distance based on fiber length is more biologically meaningful, as it reflects the direct structural basis and the cost of communication between brain regions. Longer fiber length indicates a longer physical distance for functional communication between brain regions.\nAs you hypothesized, if myelination is critically important for facilitating long-range communication, then one would expect functional coupling to be more strongly dependent on the integrity of these longer, more costly pathways. Our original analysis using geodesic distance actually served as a useful spatial measurement, but did not directly test this. To address this directly, we have conducted a new analysis as suggested.\nFirst, we incorporated DTI analysis for tractography, computed the “mean-length” of white matter fiber tract between each pair of brain regions, and defined the distance based on the group-average mean fiber tract length. We then used this axon length metric as definition of inter-regional distance for distance dependence analysis of MFC.\nSecond, to ensure robustness, we used two different neonatal brain parcellations - the UNC infant AAL atlas (82 cortical regions) and the INFANT atlas15 (211 cortical regions) to investigate the variance of MFC with spatial distance changing defined by fiber length.\nUsing the AAL atlas, for each brain region, we divided all other ipsilateral brain regions with existing white matter fiber connections into three distance subsets (short-distance, medium-distance, and long-distance) with equal counts based on non-zero group average mean-length. Detailed information regarding the fiber length for the three subsets was summarized in Table R2. The MFC was calculated for each subset independently. The line graph and brain map in Figure R12a illustrated how the mean MFC (averaged across all subjects and both hemispheres) varies across the distance subsets.\nLikewise, using the INFANT atlas (with a greater number of brain regions), we divided the ipsilateral brain regions with existing white matter fiber connections into six distance subsets (from short-distance to long-distance) with equal counts based on the non-zero group average mean-length. Detailed information regarding the fiber length for three subsets was summarized in Table R3. Then we computed the MFC for each subset, and the resulting line graph was shown in Figure R12b. The MFC shared a similar distance-dependent pattern with the results based on geodesic distance reported in our main text, with higher levels for short-distance coupling and lower levels for long-distance ones.\nOverall, the MFC pattern retained a distance-dependent profile: high coupling for short-distance connections and low coupling for long-distance connections. This characteristic is consistent across different parcellations (voxel-level: AAL atlas and INFANT atlas or surface-level: vertices) and distance definitions (cortical surface geodesic distance, or white matter fiber length between connected regions). Additionally, we found that the coupling strength of connections in the remote regions (Subset 6, defined by fiber length) can approach that of Subset 2 (Figure R12b), suggesting that, in brain network, high coupling between anatomically segregated but tightly connected core nodes supports their role in integrating whole-brain information.\nWe wish this could address your concern, and the added descriptions could further clarify the detailed method of the dominance analysis.\nFigure R12. The dependence of the MFC patterns on spatial distance defined by fiber length. a Distance dependence of MFC across three distance subsets based on the AAL atlas. b Distance dependence of MFC across six distance subsets based on the INFANT atlas.\nTable R2. The fiber length for three subsets in the AAL atlas\nThe fiber length\nSubset 1 Subset 2 Subset 3\n(mm)\nmean 42.03 59.25 67.68\nminimum 26.26 44.32 53.87\nmaximum 60.99 70.74 86.65\nTable R3. The fiber length for three subsets in the INFANT atlas\nThe fiber\nSubset 1 Subset 2 Subset 3 Subset 4 Subset 5 Subset 6\nlength (mm)\nmean 31.44 45.23 59.10 72.25 84.27 94.69\nmaximum 12.40 24.84 31.20 40.13 56.81 66.25\nminimum 62.74 71.78 87.00 97.69 100.39 100.50\nComment 5. Weak Behavioral Association: The effect size of the behavioral association is reported to be too low, which significantly weakens the clinical and biological relevance of the findings. The manuscript needs to provide a more compelling argument for the significance of this association despite the low effect size.\nResponse: We thank you for the important critique regarding the effect size of the association between neonatal MFC and 18-month behavior. We agree that the proportion of variance explained is modest, and we appreciate the opportunity to provide a more compelling interpretation of its significance within the theoretical context of our study.\nFirst, the primary aim of our work was to establish myelination covariance as a novel, microstructural axis for understanding early functional network emergence. In this framework, the behavioral association serves not as a standalone predictive model, but as a critical proof-of-concept validation. It demonstrates that this coordinated microstructural property at birth is not an epiphenomenon but carries a detectable, biologically plausible signature for later behavioral development. The statistical significance of this link, despite the temporal distance and immense intervening experiential variables, underscores its robustness as an early-life neural correlate.\nWe acknowledge that a single, neonatal brain measure cannot—and is not intended to—capture the vast complexity of behavioral development over 18 months. The modest effect size is, in fact, theoretically expected. It reflects the well-established principle that early brain organization establishes a propensity or vulnerability, which is subsequently elaborated, amplified, or mitigated by a cascade of genetic, environmental, and experiential factors. A very strong effect at this early stage would be biologically surprising, as it would leave little room for the profound developmental plasticity that characterizes infancy.\nMinor Comment 1. R-squared Correction: It is not clear from the manuscript if any correction was done for the R2 results, which is a standard procedure in such analyses to account for the number of predictors.\nResponse: Thank you for your reminder. In fact, the R2 in the research has already been corrected, and all R2 values in the results are labeled with “adj” to indicate they are adjusted. We have added a description of this standard procedure in the Methods (Section “Myelination-function coupling (MFC)”, line 693 – line 699 in manuscript):\n“In general, the coefficient of determination (R2) tends to increase when the number of parameters included in the model increases. One strategy to account for this addition is to adjust the R2 measure based on the number of predictors. Specifically, we calculated R2 for each brain region or vertex:\nadj\n𝑁−1\n𝑅2 = 1− ∙(1−𝑅2)\n𝑎𝑑𝑗 𝑁−𝑝−1\nWhere 𝑁 is the number of samples and 𝑝 is the number of predictors in the model, and is equal to 𝑝 = 2 for the dual-axis framework (sMC and gMC).”\nReferences:\n1. Honey, C. J. et al. Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A 106, 2035–2040 (2009).\n2. Horn, A., Ostwald, D., Reisert, M. & Blankenburg, F. The structural-functional connectome and the default mode network of the human brain. Neuroimage 102, 142–151 (2014).\n3. Bobba, P. S. et al. Age-related topographic map of magnetic resonance diffusion metrics in neonatal brains. Hum Brain Mapp 43, 4326–4334 (2022).\n4. Grotheer, M. et al. White matter myelination during early infancy is linked to spatial gradients and myelin content at birth. Nat Commun 13, (2022).\n5. Tau, G. Z. & Peterson, B. S. Normal development of brain circuits. Neuropsychopharmacology 35, 147–168 (2010).\n6. Sun, L. et al. Human lifespan changes in the brain’s functional connectome. Nat Neurosci 28, 891–901 (2025).\n7. Goulas, A., Betzel, R. F. & Hilgetag, C. C. Spatiotemporal ontogeny of brain wiring. Sci Adv 5, (2019).\n8. Barbas, H. & Hilgetag, C. C. From Circuit Principles to Human Psychiatric Disorders. Biol Psychiatry 93, 388–390 (2023).\n9. Jiang, L. et al. Toward neurobiological characterization of functional homogeneity in the human cortex: regional variation, morphological association and functional covariance network organization. Brain Struct Funct 220, 2485–2507 (2015).\n10. Deoni, S. C. L., Dean, D. C., Remer, J., Dirks, H. & O’Muircheartaigh, J. Cortical maturation and myelination in healthy toddlers and young children. Neuroimage 115, 147–161 (2015).\n11. Glasser, M. F., Goyal, M. S., Preuss, T. M., Raichle, M. E. & Van Essen, D. C. Trends and properties of human cerebral cortex: Correlations with cortical myelin content. Neuroimage 93, 165–175 (2014).\n12. Corrigan, N. M. et al. Myelin development in cerebral gray and white matter during adolescence and late childhood. Neuroimage 227, (2021).\n13. Seguin, C., Sporns, O. & Zalesky, A. Brain network communication: concepts, models and applications. Nat Rev Neurosci 24, 557–574 (2023).\n14. Zamani Esfahlani, F., Faskowitz, J., Slack, J., Mišić, B. & Betzel, R. F. Local structure-function relationships in human brain networks across the lifespan. Nat Commun 13, (2022).\n15. Shi, F., Salzwedel, A. P., Lin, W., Gilmore, J. H. & Gao, W. Functional Brain Parcellations of the Infant Brain and the Associated Developmental Trends. Cereb Cortex 28, 1358–1368 (2018).\n16. Kinney, H. C., Brody, B. A., Kloman, A. S. & Gilles, F. H. Sequence of central nervous system myelination in human infancy. II. Patterns of myelination in autopsied infants. J Neuropathol Exp Neurol 47, 217–234 (1988).\n17. Deoni, S. C. L., Dean, D. C., O’Muircheartaigh, J., Dirks, H. & Jerskey, B. A. Investigating white matter development in infancy and early childhood using myelin water faction and relaxation time mapping. Neuroimage 63, 1038–1053 (2012).\n18. Gilles, F. H. ., A. Leviton, . & E. C. Dooling, . The Developing Human Brain. (2013).\n19. Alexander-Bloch, A., Giedd, J. N. & Bullmore, E. Imaging structural co-variance between human brain regions. Nat Rev Neurosci 14, 322–336 (2013).\n20. Sebenius, I. et al. Structural MRI of brain similarity networks. Nat Rev Neurosci 26, 42–59 (2025).\n21. Lerch, J. P. et al. Mapping anatomical correlations across cerebral cortex (MACACC) using cortical thickness from MRI. Neuroimage 31, 993–1003 (2006).\n22. Stauffer, E. M. et al. The genetic relationships between brain structure and schizophrenia. Nat Commun 14, (2023).\n23. Váša, F. et al. Adolescent Tuning of Association Cortex in Human Structural Brain Networks. Cereb Cortex 28, 281–294 (2018).\n24. Burt, J. B. et al. Hierarchy of transcriptomic specialization across human cortex captured by structural neuroimaging topography. Nat Neurosci 21, 1251–1259 (2018).\n25. Seidlitz, J. et al. Morphometric Similarity Networks Detect Microscale Cortical Organization and Predict Inter-Individual Cognitive Variation. Neuron 97, 231-247.e7 (2018).\n26. Paquola, C. et al. Shifts in myeloarchitecture characterise adolescent development of cortical gradients. Elife 8, (2019).\n27. Sebenius, I. et al. Robust estimation of cortical similarity networks from brain MRI. Nat Neurosci 26, 1461–1471 (2023).\n28. Zhang, T. et al. Neural Plasticity Induced by Working Memory Training: Insights From Cortical Microstructure and Transcriptional Profiles. CNS Neurosci Ther 31, (2025).\n29. Wang, X. et al. Topographic, cognitive, and neurobiological profiling of the interdependent structural and functional modules of the brain. Sci Bull (Beijing) 70, 2416–2420 (2025).\n30. Zhou, Y. et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 10, (2019).\n31. Zagorski, M. et al. Decoding of position in the developing neural tube from antiparallel morphogen gradients. Science 356, 1379–1383 (2017).\n32. Briscoe, J. & Ericson, J. Specification of neuronal fates in the ventral neural tube. Curr Opin Neurobiol 11, 43–49 (2001).\n33. Karzbrun, E. et al. Human neural tube morphogenesis in vitro by geometric constraints. Nature 2021 599:7884 599, 268–272 (2021).\n34. Arai, K. & Lo, E. H. Gliogenesis. Primer on Cerebrovascular Diseases: Second Edition 91–95 (2017) doi:10.1016/B978-0-12-803058-5.00018-7.\n35. Sohur, U. S., Emsley, J. G., Mitchell, B. D. & Macklis, J. D. Neural Stem Cells. Encyclopedia of Stress 832–840 (2007) doi:10.1016/B978-012373947-6.00574-2.\n36. Vieira, M. S. et al. Decoding epigenetic cell signaling in neuronal differentiation. Semin Cell Dev Biol 95, 12–24 (2019).\n37. Gelfo, F., Mandolesi, L., Serra, L., Sorrentino, G. & Caltagirone, C. The Neuroprotective Effects of Experience on Cognitive Functions: Evidence from Animal Studies on the Neurobiological Bases of Brain Reserve. Neuroscience 370, 218–235 (2018).\n38. Romero, D. M., Bahi-Buisson, N. & Francis, F. Genetics and mechanisms leading to human cortical malformations. Semin Cell Dev Biol 76, 33–75 (2018).\n39. Camacho-Arroyo, I., Piña-Medina, A. G., Bello-Alvarez, C. & Zamora-Sánchez, C. J. Sex hormones and proteins involved in brain plasticity. Vitam Horm 114, 145–165 (2020).\n40. Power, J. D. et al. Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage 84, 320–341 (2014).\n41. Fitzgibbon, S. P. et al. The developing Human Connectome Project (dHCP) automated resting-state functional processing framework for newborn infants. Neuroimage 223, (2020).\n42. Andersson, J. L. R. et al. Towards a comprehensive framework for movement and distortion correction of diffusion MR images: Within volume movement. Neuroimage 152, 450–466 (2017).\n43. Ciric, R. et al. Benchmarking of participant-level confound regression strategies for the control of motion artifact in studies of functional connectivity. Neuroimage 154, 174–187 (2017).\n44. Kim, C. N., Shin, D., Wang, A. & Nowakowski, T. J. Spatiotemporal molecular dynamics of the developing human thalamus. Science 382, (2023).\n45. Alhesain, M. et al. Development of the early fetal human thalamus: from a protomap to emergent thalamic nuclei. Frontiers in Neuroanatomy 19, (2025).\n46. Zhang, J. et al. Cortical and subcortical mapping of the human allostatic-interoceptive system using 7 Tesla fMRI. Nat Neurosci 28, 2380–2391 (2025).\n47. Hayashi, S. et al. Maturation of Complex Synaptic Connections of Layer 5 Cortical Axons in the Posterior Thalamic Nucleus Requires SNAP25. Cereb Cortex 31, 2625–2638 (2021).\n48. Ruiz-Cabrera, S. et al. Expansion modes of primate nervous system structures in the light of the Prosomeric Model. Frontiers in Mammal Science 2, 1241573 (2023).\nDear Editor and Reviewers,\nWe sincerely thank the editor and the reviewers. During this revision, we have carefully considered the constructive comments raised by the Reviewer 2 as well as all the editorial requests.\nFor convenience in checking the changes made in the revised manuscript, we have marked these changes with ‘red’ color. We wish our revised manuscript would be more readable than the previous version and reach the high standard requirement of this journal.\nThank you once again for your time and expertise.\nBest wishes,\nYours sincerely,\nWeijin Liu & Yuanyuan Chen\nOn behalf of all authors\nREVIEWER COMMENTS AND AUTHOR RESPONSE\nReviewer #1 (Remarks to the Author):\nThe authors have adequately addressed my comments and concerns.\nResponse: Thank you for your supportive comment.\nReviewer #2 (Remarks to the Author):\nI would like to thank the authors for providing a solid revision that carefully addresses many of my concerns. While most issues have been resolved, I still have two additional major comments and one minor comment.\nComment 1. Throughout the manuscript, the authors continue to ascribe a causal relationship between myelination and functional coupling (e.g., “offering mechanistic insights into how coordinated myelination may drive the initial assembly of functional pathways”). From my understanding, the presented data are correlational in nature. An opposite causal relationship—where coordinated functional activity shapes myelination—is equally plausible, particularly given the well-established activity dependence of myelination. I therefore strongly recommend revising the manuscript to avoid causal language and to frame these findings in explicitly correlational terms.\nResponse: We thank the reviewer for this essential correction and fully agree with the point. Our intention was to propose a plausible mechanistic hypothesis which could supplement white matter only as a structural foundation for early functional connectivity development, but as what you mentioned that our language sometimes overstepped the correlational nature of the evidence. We have now thoroughly revised the manuscript to eliminate causal language and consistently adopt a correlational framework.\nFirstly, in Abstract & Discussion, we have replaced terms like ”drive”, and “scaffold” with more precise language such as “covaries with”, “couples with”, and “may reflect.” For instance, the phrase quoted by the reviewer has been revised to: “offering mechanistic hypotheses for how coordinated myelination covaries with the initial assembly of functional pathways.”\nSecondly, we carefully concluded and interpreted our findings as association or a framework to describe a coordinated, parallel development or co-variation between intracortical myelination covariance and functional synchronization, rather than a unidirectional driving relationship. We explicitly stated in the Discussion (line 610 - 615) that the observed associations do not establish causality and that activity-dependent myelination remains a highly plausible, reciprocal mechanism.\n“…while associations between myelination patterns and functional connectivity are demonstrated, the causal mechanisms linking myelination covariance to functional network emergence remain unclear, suggesting that other neurodevelopmental processes may also contribute. Experimental studies, including animal models, will be crucial to better understand the causality of these relationships.”\nWe believe these revisions bring the manuscript into full alignment with the evidence and the reviewer’s valid critique.\nComment 2. I appreciated the additional analysis of structural connectivity and found it very interesting that structural connectivity explains additional variance in functional connectivity. Even though this effect is relatively small, it adds important depth to the developmental relationships proposed by the authors. For this reason, I feel it is inappropriate to relegate this analysis to the supplemental material and would instead suggest devoting an additional Results subsection to these findings. Moreover, the authors’ interpretation that “These findings suggest that during infancy, the integration of brain function is not initially dominated by long-range white matter connections but is grounded in the synchrony of intracortical microstructure that reflects shared developmental trajectories” does not fully reflect the reported results, as both white- and gray-matter structural measures relate to functional connectivity. Finally, I could not find sufficient methodological details describing how the diffusion MRI analyses were performed. Please ensure that the Methods section fully describes all analyses, including those currently presented only in the Supplementary Materials.\nResponse: We thank the reviewer for highlighting the importance of these analyses and for the helpful suggestions to improve their presentation and interpretation.\nFirstly, we agree that the comparison between gray matter MFC and white matter structure-function coupling (SFC) is a core part of the story. As suggested, we have now moved the results about white matter based SFC analysis from the Supplementary materials to the Results as a new section titled “Limited contribution from the white matter connectivity for the coupling.” :\n“Having established that structural covariance of gray matter myelination is associated with the emergence of functional connectivity, it remained possible that white matter fiber connections also play a role therein. Therefore, we incorporated DTI data for fiber tract reconstruction, calculated the fiber tract connectivity density between brain regions to construct the structural connectivity network (SCN) based on the UNC infant AAL atlas. Subsequently, the fully weighted SCN was used to model the coupling relationship between white matter tracts and functional connectivity. The coupling was calculated using the same method as MFC, adopting a linear regression model. Averaged map of atlas-level individual MFC and SFC were illustrated (Supplementary Fig. 5a and b), showing that a moderate level of correlation between them (corr = 0.48, p < 0.0001; Supplementary Fig.5c). MFC map showed much higher level than SFC across the whole brain (the mean MFC is 0.22, and the mean SFC is 0.09), except for the cuneate cortex, with the peak increment up to 9 times (Supplementary Fig. 5d, the coupling increment map of MFC over SFC). The dominant regions of SFC were only distributed in the medial cortex, including the primary visual cortex and primary motor cortex. Furthermore, we introduced structural connectivity (SC) as an additional predictive variable into the MFC dual-axis framework. For the majority of brain regions, the introduction of SC profile barely improved the prediction strength for functional connectivity, particularly in the dorsolateral prefrontal cortex (Supplementary Fig. 5e, MFC map with SC included). Regions with an increment exceeding 20% were concentrated in the posterior cingulate cortex, calcarine cortex, and cuneus (Supplementary Fig. 5f, Coupling increment map of SC-included MFC over MFC without SC variable).”\nSecondly, for the sentence in abstract, we have rephased it as:\n“…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, ….”\nThirdly, we apologize for the omission of methodological detail in the main text. We have now moved the relevant method descriptions from Supplementary materials to the Methods as a new section “White matter based structural connectivity network”:\n“Structural connectivity (SC) between each pair of regions was calculated as the fiber tract connectivity density using DTI imaging data for fiber tract reconstruction based on the ROI-level atlas. We utilized two alternative volume-based infant atlas: the UNC infant AAL atlas91, comprising 8 subcortical structures and 82 cortical structures. The reconstruction of the diffusion data was performed in native space with DSI-studio (https://dsi-studio.labsolver.org/). After the generalized q-sampling imaging (GQI) with a diffusion sampling length ratio of 1.25, the whole-brain fiber tracking was conducted with the parameters including the angular cutoff of 60 degrees, step size of 1.0 mm, minimum length of 30 mm, and maximum length of 300 mm. The whole-brain fiber tracking process was performed with the FACT algorithm until 1,000,000 streamlines were reconstructed for each individual. The fiber tract connectivity density was quantified as the number of streamlines normalized by volumes of paired regions, yielding a 90 × 90 (90 parcellations with 82 cortical regions and 8 subcortical regions) matrix encoded the structural connectivity network (SCN) For each participant. Then, Sparse structural connections were converted into a communicability matrix to avoid the loss of structural information, resulting in a fully weighted form of the SCN based on the communication model92,93.”\nWe are confident that these revisions, inspired by the reviewer’s insightful comments, have substantially strengthened the clarity, accuracy, and balance of our manuscript. We deeply appreciate the time and expertise invested in helping us improve this work.\nComment 3. There are additional publications on gray/white matter myelin coupling that the authors may wish to consider citing:\nhttps://doi.org/10.1007/s00429-025-03054-9\nhttps://doi.org/10.1101/2025.07.07.663449\nhttps://doi.org/10.1038/s41467-026-68714-8\nResponse: Thanks very much for your recommended publications and we have cited them at relevant context in the main text.\nReviewer #3 (Remarks to the Author):\nAll my concerns have been addressed.\nResponse: Thank you for your supportive comment."}]

Decision letter