Abstract
Functional organization of the human cerebral cortex is highly constrained by underlying brain structures, but how functional activity is associated with different brain structures during development is not clear, especially at the neonatal stage. Since long-range functional connectivity is far from mature in the dynamically developing neonatal brain, it is of great scientific significance to investigate the relationship between different structural and functional features at the local level. To this end, for the first time, correlation and regression analyses were performed to examine the relationship between cortical morphology, cortical myelination, age, and local brain functional activity, as well as functional connectivity strength using high-resolution structural and resting-state functional MRI data of 177 neonates (29–44 postmenopausal weeks, 98 male and 79 female) from both static and dynamic perspectives. We found that cortical myelination was most strongly associated with local brain functional activity across the cerebral cortex than other cortical structural features while controlling the age effect. These findings suggest the crucial role of cortical myelination in local brain functional development at birth, providing valuable insights into the fundamental biological basis of functional activity at this early developmental stage.
Introduction
How various brain functional activities arise from the complex cortical structures is a fundamental problem in neuroscience. Spatially heterogeneous associations between structural connectivity or covariance (using cortical morphology) and functional interactions in adults have been found (Honey et al. 2009; Smith et al. 2019; Vázquez-Rodríguez et al. 2019; Gu et al. 2021). However, to date, how the structure–function coupling is established and evolves during early brain development is poorly understood, especially in neonates.
The human brain develops extremely rapidly in the first postnatal months, with enlarged brain volumes (Gilmore et al. 2007; Gilmore et al. 2020) and cortical thickness (Li et al. 2015a; Li et al. 2015b), more complex gyrification (Li et al. 2014) and functional connections (Gao et al. 2015a; Keunen et al. 2017), and increased myelin content (Dubois et al. 2014). However, how functional activities arise from brain structures during this dramatically developing neonatal stage is not clear. It is of great scientific importance to explore the brain structure–function relationship at birth, since the neonatal brain is largely immature, which remarkably differs from the adult brain, and the structure–function correspondence is mainly genetically determined with minimal environmental influences at this stage. In the literature, cortical morphology has been extensively used to investigate its relationship with functional brain activity and connectivity (Biagioni et al. 2007; Geng et al. 2017), as well as cognitive ability and behavior outcomes during infancy (Shaw et al. 2006; Nam et al. 2015; Adeli et al. 2019; Girault et al. 2020). However, the biological rationale of using morphological features and whether there is a better predictor for brain functions in infants remains to be elucidated.
A systematic investigation of the potential correspondence between structure and function in the neonatal brain is highly needed, considering previous studies’ limitations (Biagioni et al. 2007; Geng et al. 2017; Huntenburg et al. 2017) and neonatal characteristics. Although resting-state functional connectivity or networks have been widely investigated in adult studies (Paquola et al. 2019; Smith et al. 2019), long-range functional connections have not well developed at birth (Gao et al. 2015b), indicating local functional features such as local brain activity are a better choice to study in neonates. Previous adult and infant studies mainly employed a single or only a few brain structural and functional features (Biagioni et al. 2007; Geng et al. 2017; Huntenburg et al. 2017). However, for the neonatal brain, it would be more promising to use multiple features to capture more aspects and complementary information of brain development (Fenchel et al. 2020). For example, myelination, an important feature reflecting the conduction velocity of nerve impulses along axons and potentially relating to cortical functions (Fornari et al. 2007), and whose maturation lasts from the third trimester to adolescents (Nave and Werner 2014), is rarely leveraged but crucial for early brain development (Deoni et al. 2011).
In this paper, for the first time, we systemically characterize the relationship between cortical structural and local functional features as well as functional connectivity-based features of the neonatal brain. First, we obtain cortical structural (morphological features and myelination) and functional indices (static and dynamic local brain activity and functional connectivity strength) of 177 healthy neonates. We then investigate the relationship between age, structural and functional indices using correlation analysis and regression analysis. The study helps answer how and which structural features are closely associated with functional brain activity in neonates and sheds new light on brain structure–function correspondence during infancy.
Materials and methods
Subjects
Subjects included in this study were from the second release of an international publicly available database: developing human connectome project (Hughes et al. 2017) (dHCP, https://data.developingconnectome.org). A total of 177 healthy subjects were picked out from 505 neonates according to inclusion criteria: (i) T1w, T2w, and resting-state functional images are in place at the same time; (ii) cerebral cortex is covered in MR images and there are no image quality problems by visual inspection; (iii) the extracted functional period meets the criteria of head motion. The 177 subjects included 51 preterm infants (post-menstrual age < 37 weeks at birth) and 126 full-term neonates and consisted of 98 males and 79 females. These subjects were scanned during natural sleep without sedation in postmenopausal 29–44 weeks, and the mean scan age was 39.54(±3.03) weeks (details in Fig. 1). The gap between scan and birth was 1.94(±3.06) weeks, and 105 subjects were scanned within one week after birth. There were 31 preterms at scan. Preterm infants were included in the current study to examine the structure–function correspondence in a larger age span and to help reduce the bias caused by fewer sampling points at a younger age in the process of model estimation. A validation analysis was performed without preterm infants (at scan) to examine the reliability of the relationship between structural and functional indices in a narrower age range.
Birth age and scan age of subjects in this study.
Data acquisition
Neonatal MR images were acquired when neonates were sleeping using Philips Achieva 3 T scanner and a neonatal 32 channel phased array head coil at the Evelina Neonatal Imaging Centre, St Thomas’ Hospital, London (Hughes et al. 2017). Multiband 9× accelerated echo-planar imaging with 2.15-mm isotropic resolution (repeated time (TR) = 392 ms, time echo (TE) = 38 ms) was used to obtain functional images for 15 min. Multi-slice fast spin-echo structural images (in-plane resolution = 0.8 × 0.8 mm2, 1.6-mm slices overlapped by 0.8 mm) were acquired with the following parameters: for T1 weighted images, TR = 4,795 ms, TE = 8.7 ms; for T2-weighted images, TR = 12,000 ms, TE = 156 ms.
Image processing
We used the structural and functional MR images preprocessed by the dHCP pipeline. Preprocessing steps of structural images included bias correction, skull-stripping, tissue segmentation, T1w-T2w alignment, and surface extraction (Makropoulos et al. 2018). Preprocessing steps of functional images comprised the correction of susceptibility distortions, motion correction, structural-functional alignment, high-pass filter (150 s), and ICA denoising (Fitzgibbon et al. 2020). Band-filtering (0.01–0.1 Hz) was applied to functional images using DPABI (Yan et al. 2016). Global signal regression was performed on blood-oxygen-level-dependent (BOLD) signals using RESTplus (Jia et al. 2019) as a validation analysis.
The reconstructed cortical surfaces from dHCP-released packages were mapped onto a sphere to achieve a spherical representation of the surface for inter-subject registration. Since the dHCP pipeline used a multidimensional scaling method to achieve the spherical mapping, the generated spherical surfaces sit in a unique space that may not be consistent with most existing atlases. To be consistent with the atlas space, we used our own post-processing steps. Specifically, we first mapped the white surface to a sphere and then aligned the generated spherical surfaces onto the neonatal dedicated cortical surface atlases (Wu et al. 2019) using the FreeSurfer toolkit (http://surfer.nmr.mgh.harvard.edu/). We choose this atlas because it provided the week-wise cortical surface atlases from 39 to 44 gestational weeks. Considering the rapid brain development during the neonatal stage, a high-definition atlas in time can better chart the cerebral cortex with clear folding patterns and lead to better registration accuracies. After the registration, we linearly resampled the registered spherical surface to be 10,242 vertices for each hemisphere with a uniform sphere tessellation.
Functional images were extracted given that big head movement occurs occasionally during infants’ natural sleep (>5 mm). The criteria of extracted periods were as follows: (i) continuous periods; (ii) framewise displacement (FD) smaller than 1 mm in the period; (iii) at least 3 min, i.e. 459 volumes and (iv) the minimal mean FD in all possible periods to be extracted. We did not perform data censoring or concatenation to preserve the temporal continuity of BOLD signals, which may influence the following sliding-window analyses (Hutchison et al. 2013). The use of at least 3-min functional periods refers to previous studies (Wen et al. 2019), and 459 volumes are sufficient when compared to other infant studies (Gao et al. 2015a) or dynamic functional connectivity studies (Allen et al. 2014). Finally, functional periods from 177 subjects were used in this study with an average of 484 volumes (3.16 min), and the mean FD was 0.12(±0.03) mm, which was smaller compared to other infant studies (Yin et al. 2020).
Structural indices
Structural indices were made up of cortical morphology and cortical myelination (CMY). Cortical morphology was often described using cortical thickness (CT), surface area (SA), and local gyrification index (LGI). CT is defined as the distance between inner and outer cortical surfaces (Fischl and Dale 2000), while SA is defined as the area of a vertex or a region on the surface mesh. LGI measures the complexity of cortical folding, which creates a sphere with a 15-mm radius based on each vertex and then calculates the ratio between the surface area within the sphere and the area of the great circle of the sphere (Schaer et al. 2008). CMY measures the myelination of axons in the cerebral cortex, reflecting the efficiency of signal transmission between neurons, and it is defined as the cortical part of the myelin content map (T1w/T2w ratio map (Glasser and van Essen 2011)). The mid-thickness surface is the middle surface between inner and outer surfaces, which can be obtained by calculating the mean Euclidean distance between the corresponding vertices of inner and outer surfaces. CT and myelin content maps were provided by the dHCP pipeline, CT was resampled, and the myelin map was mapped to the mid-thickness surface to obtain CMY. To keep consistency with functional indices, which also relied on the method of volume-to-surface mapping using the software Workbench (https://www.humanconnectome.org/software/connectome-workbench), the myelin map and mid-thickness surface were transformed from T2w to functional space. SA and LGI were calculated based on the mid-thickness surface, and SA was obtained using Workbench. Gaussian smooth with an 8-mm Full-width Half Maximum (FWHM) kernel on the mid-thickness surface was performed after the structural indices were obtained. To reduce the impact of arbitrary choice of kernel size, we additionally used a smaller size (4 mm, the same as in dHCP pipeline) to perform the analysis at the whole brain level.
Functional indices
We used various functional brain measures aiming to capture different aspects of the neonatal brain: static and dynamic local-brain and connectivity-based features. Amplitude of low-frequency fluctuations (ALFF) (Zang et al. 2007) and fractional ALFF (fALFF)(Zou et al. 2008) were leveraged to describe local brain functional activity in the frequency domain. Regional homogeneity (ReHo) and functional connectivity strength (FCS) were used to measure local connectivity and global connectivity respectively. Corresponding dynamic indices (d-ALFF, d-fALFF, d-ReHo, and stdFCS) were also obtained here to explore the temporal variability of the indices during the scan with the technique of sliding window.
Specifically, ALFF measures the amplitude fluctuation of spontaneous neural activity in low frequency (Here, we used 0.01–0.1 Hz). fALFF measures the ratio of spontaneous neural activity in low frequency to the entire frequency (0.0067–2.55 Hz). ReHo measures the regional homogeneity of spontaneous neural activity, which is obtained by calculating Kendall’s concordance coefficient of timeseries between one voxel and its nearest 26 neighbors. FCS measures the average strength of functional connectivity of a given region with the rest of the brain. The Fisher’s r-to-z transformation was applied to increase the normality of the distribution of correlation coefficients. Connectivity thresholds were set to exclude spurious connections, and absolute values were used to obtain the strength of both positive and negative connectivity above the threshold. For each hemisphere, the threshold was determined using a one-sample T-test on the functional connectivity profile for each region, the significance was corrected by family-wise error (FWE) correction (|$P<0.05$|), and the thresholds were averaged across regions. For the analysis based on vertex, the threshold |${r}_{th}$| was 0.1990 for the left hemisphere and 0.1989 for the right; for the analysis based on parcellation, |${r}_{th}$| was 0.1547 for the left hemisphere and 0.1566 for the right.
The computation of dynamic indices combined the static indices mentioned above and the sliding window technique, which contained: (i) split the timeseries of each voxel into a series of time windows with 40 s as the window length and 1 TR (0.392 s) as the step size; (ii) calculate ALFF, fALFF, ReHo, and FCS in each time window; (iii) calculate the mean value and standard deviation of indices across time windows; (iv) For ALFF, fALFF, and ReHo, dividing the standard deviation by the mean (i.e. coefficient of variation) across time windows is defined as the corresponding dynamic indices (d-ALFF, d-fALFF, and d-ReHo); For FCS, the standard deviation across windows is defined as dynamic FCS (stdFCS). The choice of 40 s as window length followed suggestions by (Hutchison et al. 2013), who proposed 30–60 s as proper lengths. To reduce the impact of arbitrary choice of window length, larger and smaller windows (i.e. 50 and 30 s) were considered in the validation analysis. The choice of 1TR step was consistent with numerous previous studies (Allen et al. 2014). There were 358–664 time windows and 383 on average due to different numbers of time points among subjects.
Functional indices except static and dynamic FCS were calculated by DPASF5.0 (Yan et al. 2016). Timeseries and local features of functional images in volume space were mapped onto the mid-thickness surface using Workbench and then smoothed using an 8-mm FWHM kernel as was done on structural indices.
Multiscale analysis
This study was based on the whole brain, vertex, and parcellation. For each index, the whole-brain level of each subject was the average between the left and right hemispheres except SA, which was defined as the sum of the two hemispheres. We applied the method of functional parcellation that we proposed in our previous work (Zhao et al. 2020) on dHCP MRI data to acquire 190 cortical regions in the left hemisphere and 193 in the right since there was no proper surface-based neonate-specific functional atlas. For CT, LGI, CMY, static and dynamic ALFF, fALFF, and ReHo, the regional indices were the average value of vertices in the same region. For SA, the regional SA was the sum of SA of vertices in the same region. For FCS and stdFCS, average the timeseries of vertices in the same region to acquire the regional timeseries, and then calculate FCS and stdFCS above the connectivity threshold |${r}_{th}.$| The calculation of regional |${r}_{th}$| was the same as that of vertex-level.
Statistical analysis
Correlation analysis
The correlations between age/structural indices and functional indices were explored from the whole-brain and vertex levels. The Pearson correlation coefficients between age/structural indices and functional indices were calculated, and the corresponding |$p$| values were corrected using Bonferroni correction. For the whole-brain level, there were 78 (13 × 12/2 = 78) times comparisons, which was the number of elements in the lower triangular correlation matrix of 13 indices. For the vertex level, there were around 9,100 (the number of cortical vertices in the left or right hemisphere) times comparisons, and cluster-wise FDR correction was applied with a minimum area of 30 mm2 in a cluster using FreeSurfer. The significance level was 0.05 for the whole study.
Regression analysis
The study of the relationship between structural and functional development virtually examined the relationship between two kinds of equations: |$struct={\beta}_s age+{\varepsilon}_s$| and |$func={\beta}_f age+{\varepsilon}_f$|, where |$age$| indicates scan age, |$struct$| indicates structural indices, and |$func$| indicates functional indices. And the relationship between these equations could be mathematically transformed as |$func={\beta}_1 age+{\beta}_2 struct+\varepsilon .$| Since the neonatal brain is developing rapidly, all the structural and functional indices are correlated with age. Including age in the regression model as an independent variable can also control the impact of common age dependence when exploring the relationship between structural and functional indices.
Three types of multiple regression models were estimated to fit the relationship between structural and functional development:
- (1) Linear model:$$ f=\alpha +{\beta}_1 age+{\beta}_2 CT+{\beta}_3 SA+{\beta}_4 LGI+{\beta}_5 CMY+\varepsilon $$
- (2) Logarithmic model:$$ f=\alpha +{\beta}_1 age+{\beta}_2\log (age)+{\beta}_3 CT+{\beta}_4 SA+{\beta}_5 LGI+{\beta}_6 CMY+\varepsilon $$
- (3) Quadratic model:$$ f=\alpha +{\beta}_1 age+{\beta}_2{age}^2+{\beta}_3 CT+{\beta}_4 SA+{\beta}_5 LGI+{\beta}_6 CMY+\varepsilon $$
where |$f$| represents one functional index each time, and |$age$| represents the scan age. Post-menstrual age at birth (BirthAge), birth weight (weight), and gender were included in models as covariates at the global level, and only the significant covariates (|$P<0.05$|) were included in further studies. Since residual motion-related artifacts were separated from BOLD signals in the process of ICA denoising by the dHCP pipeline (Fitzgibbon et al. 2020), we did not include motion-related variables in the regression models. Global signal regression was applied to BOLD signals, and then the same regression analyses at the global and vertex level were performed as a validation analysis.
Adjusted R square (|$adj{R}^2$|) was used to assess the goodness-of-fit of the models and the model with the biggest |$adj{R}^2$| (average |$adj{R}^2$| for the vertex and region levels) was the best model. At the whole-brain level, backward elimination stepwise regression was also performed for each functional index, which included all predictors at first and eliminated the insignificant predictors (|$P<0.05$|), and the final models were called partial models. Cohen’s f2 values were also reported to evaluate the effect size of the models.
Multiple comparison corrections were performed on the regression results in vertex and region levels by cluster-wise FDR correction with a minimum area of 30 mm2 for the former using FreeSurfer and FDR correction for the latter using MATLAB (both |$P<0.05$|), and the number of comparisons was the number of cortical vertices or parcellations in each hemisphere.
Results
Correlations between age, structural, and functional indices
A variety of structural (CT, SA, LGI, and CMY) and functional indices (static and dynamic ALFF, fALFF, ReHo, and FCS) were used to explore the relationship between brain structure and function of neonates. The spatial distribution maps of indices were shown in a Supplementary Fig. S1. Pearson correlation coefficients were computed between scan age (age), structural and functional indices for the whole-brain and vertex levels. Age and CMY showed the highest correlation with other features (mean r = 0.55), while CT showed the lowest (mean r = 0.19)(as shown in Table 1 and Supplementary Fig. S2). Among functional indices, d-ReHo was the only one negatively correlated with age or structural indices. Correlations were all significant after Bonferroni correction (maximal corrected |$P$| = 0.0420) except for those between CT/d-ReHo and certain indices.
Pearson correlation coefficients between age, structural, and functional indices for the whole brain.
| Age | CT | SA | LGI | CMY | ALFF | fALFF | ReHo | d-ALFF | d-fALFF | d-ReHo | FCS | stdFCS | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Age | 1.00 | – | – | – | – | – | – | – | – | – | – | – | – |
| CT | 0.34*** | 1.00 | – | – | – | – | – | – | – | – | – | – | – |
| SA | 0.92*** | 0.32** | 1.00 | – | – | – | – | – | – | – | – | – | – |
| LGI | 0.94*** | 0.13 | 0.88*** | 1.00 | – | – | – | – | – | – | – | – | – |
| CMY | 0.84*** | 0.24 | 0.77*** | 0.85*** | 1.00 | – | – | – | – | – | – | – | – |
| ALFF | 0.50*** | 0.25 | 0.42*** | 0.44*** | 0.49*** | 1.00 | – | – | – | – | – | – | – |
| fALFF | 0.57*** | 0.09 | 0.53*** | 0.59*** | 0.66*** | 0.59*** | 1.00 | – | – | – | – | – | – |
| ReHo | 0.59*** | 0.08 | 0.54*** | 0.63*** | 0.67*** | 0.50*** | 0.95*** | 1.00 | – | – | – | – | – |
| d-ALFF | 0.50*** | 0.25 | 0.48*** | 0.46*** | 0.53*** | 0.38*** | 0.67*** | 0.64*** | 1.00 | – | – | – | – |
| d-fALFF | 0.46*** | 0.32** | 0.46*** | 0.40*** | 0.45*** | 0.29*** | 0.48*** | 0.45*** | 0.96*** | 1.00 | – | – | – |
| d-ReHo | −0.26* | 0.14 | −0.21 | −0.31** | −0.28* | −0.19 | −0.51*** | −0.56*** | 0.15 | 0.35*** | 1.00 | – | – |
| FCS | 0.32*** | −0.09 | 0.31** | 0.38*** | 0.37*** | 0.27* | 0.71*** | 0.74*** | 0.53*** | 0.38*** | −0.45*** | 1.00 | – |
| stdFCS | 0.34*** | −0.08 | 0.33*** | 0.41*** | 0.40*** | 0.27* | 0.72*** | 0.76*** | 0.65*** | 0.50*** | −0.31** | 0.95*** | 1.00 |
| Mean | 0.55 | 0.19 | 0.51 | 0.54 | 0.55 | 0.38 | 0.59 | 0.59 | 0.52 | 0.46 | 0.31 | 0.46 | 0.48 |
| Age | CT | SA | LGI | CMY | ALFF | fALFF | ReHo | d-ALFF | d-fALFF | d-ReHo | FCS | stdFCS | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Age | 1.00 | – | – | – | – | – | – | – | – | – | – | – | – |
| CT | 0.34*** | 1.00 | – | – | – | – | – | – | – | – | – | – | – |
| SA | 0.92*** | 0.32** | 1.00 | – | – | – | – | – | – | – | – | – | – |
| LGI | 0.94*** | 0.13 | 0.88*** | 1.00 | – | – | – | – | – | – | – | – | – |
| CMY | 0.84*** | 0.24 | 0.77*** | 0.85*** | 1.00 | – | – | – | – | – | – | – | – |
| ALFF | 0.50*** | 0.25 | 0.42*** | 0.44*** | 0.49*** | 1.00 | – | – | – | – | – | – | – |
| fALFF | 0.57*** | 0.09 | 0.53*** | 0.59*** | 0.66*** | 0.59*** | 1.00 | – | – | – | – | – | – |
| ReHo | 0.59*** | 0.08 | 0.54*** | 0.63*** | 0.67*** | 0.50*** | 0.95*** | 1.00 | – | – | – | – | – |
| d-ALFF | 0.50*** | 0.25 | 0.48*** | 0.46*** | 0.53*** | 0.38*** | 0.67*** | 0.64*** | 1.00 | – | – | – | – |
| d-fALFF | 0.46*** | 0.32** | 0.46*** | 0.40*** | 0.45*** | 0.29*** | 0.48*** | 0.45*** | 0.96*** | 1.00 | – | – | – |
| d-ReHo | −0.26* | 0.14 | −0.21 | −0.31** | −0.28* | −0.19 | −0.51*** | −0.56*** | 0.15 | 0.35*** | 1.00 | – | – |
| FCS | 0.32*** | −0.09 | 0.31** | 0.38*** | 0.37*** | 0.27* | 0.71*** | 0.74*** | 0.53*** | 0.38*** | −0.45*** | 1.00 | – |
| stdFCS | 0.34*** | −0.08 | 0.33*** | 0.41*** | 0.40*** | 0.27* | 0.72*** | 0.76*** | 0.65*** | 0.50*** | −0.31** | 0.95*** | 1.00 |
| Mean | 0.55 | 0.19 | 0.51 | 0.54 | 0.55 | 0.38 | 0.59 | 0.59 | 0.52 | 0.46 | 0.31 | 0.46 | 0.48 |
“mean” represents the average value of the absolute correlation coefficient of one index with others. *: corrected |$P<0.05$|, **: corrected |$P<0.01$|, ***: corrected |$P<0.001$|.
Pearson correlation coefficients between age, structural, and functional indices for the whole brain.
| Age | CT | SA | LGI | CMY | ALFF | fALFF | ReHo | d-ALFF | d-fALFF | d-ReHo | FCS | stdFCS | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Age | 1.00 | – | – | – | – | – | – | – | – | – | – | – | – |
| CT | 0.34*** | 1.00 | – | – | – | – | – | – | – | – | – | – | – |
| SA | 0.92*** | 0.32** | 1.00 | – | – | – | – | – | – | – | – | – | – |
| LGI | 0.94*** | 0.13 | 0.88*** | 1.00 | – | – | – | – | – | – | – | – | – |
| CMY | 0.84*** | 0.24 | 0.77*** | 0.85*** | 1.00 | – | – | – | – | – | – | – | – |
| ALFF | 0.50*** | 0.25 | 0.42*** | 0.44*** | 0.49*** | 1.00 | – | – | – | – | – | – | – |
| fALFF | 0.57*** | 0.09 | 0.53*** | 0.59*** | 0.66*** | 0.59*** | 1.00 | – | – | – | – | – | – |
| ReHo | 0.59*** | 0.08 | 0.54*** | 0.63*** | 0.67*** | 0.50*** | 0.95*** | 1.00 | – | – | – | – | – |
| d-ALFF | 0.50*** | 0.25 | 0.48*** | 0.46*** | 0.53*** | 0.38*** | 0.67*** | 0.64*** | 1.00 | – | – | – | – |
| d-fALFF | 0.46*** | 0.32** | 0.46*** | 0.40*** | 0.45*** | 0.29*** | 0.48*** | 0.45*** | 0.96*** | 1.00 | – | – | – |
| d-ReHo | −0.26* | 0.14 | −0.21 | −0.31** | −0.28* | −0.19 | −0.51*** | −0.56*** | 0.15 | 0.35*** | 1.00 | – | – |
| FCS | 0.32*** | −0.09 | 0.31** | 0.38*** | 0.37*** | 0.27* | 0.71*** | 0.74*** | 0.53*** | 0.38*** | −0.45*** | 1.00 | – |
| stdFCS | 0.34*** | −0.08 | 0.33*** | 0.41*** | 0.40*** | 0.27* | 0.72*** | 0.76*** | 0.65*** | 0.50*** | −0.31** | 0.95*** | 1.00 |
| Mean | 0.55 | 0.19 | 0.51 | 0.54 | 0.55 | 0.38 | 0.59 | 0.59 | 0.52 | 0.46 | 0.31 | 0.46 | 0.48 |
| Age | CT | SA | LGI | CMY | ALFF | fALFF | ReHo | d-ALFF | d-fALFF | d-ReHo | FCS | stdFCS | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Age | 1.00 | – | – | – | – | – | – | – | – | – | – | – | – |
| CT | 0.34*** | 1.00 | – | – | – | – | – | – | – | – | – | – | – |
| SA | 0.92*** | 0.32** | 1.00 | – | – | – | – | – | – | – | – | – | – |
| LGI | 0.94*** | 0.13 | 0.88*** | 1.00 | – | – | – | – | – | – | – | – | – |
| CMY | 0.84*** | 0.24 | 0.77*** | 0.85*** | 1.00 | – | – | – | – | – | – | – | – |
| ALFF | 0.50*** | 0.25 | 0.42*** | 0.44*** | 0.49*** | 1.00 | – | – | – | – | – | – | – |
| fALFF | 0.57*** | 0.09 | 0.53*** | 0.59*** | 0.66*** | 0.59*** | 1.00 | – | – | – | – | – | – |
| ReHo | 0.59*** | 0.08 | 0.54*** | 0.63*** | 0.67*** | 0.50*** | 0.95*** | 1.00 | – | – | – | – | – |
| d-ALFF | 0.50*** | 0.25 | 0.48*** | 0.46*** | 0.53*** | 0.38*** | 0.67*** | 0.64*** | 1.00 | – | – | – | – |
| d-fALFF | 0.46*** | 0.32** | 0.46*** | 0.40*** | 0.45*** | 0.29*** | 0.48*** | 0.45*** | 0.96*** | 1.00 | – | – | – |
| d-ReHo | −0.26* | 0.14 | −0.21 | −0.31** | −0.28* | −0.19 | −0.51*** | −0.56*** | 0.15 | 0.35*** | 1.00 | – | – |
| FCS | 0.32*** | −0.09 | 0.31** | 0.38*** | 0.37*** | 0.27* | 0.71*** | 0.74*** | 0.53*** | 0.38*** | −0.45*** | 1.00 | – |
| stdFCS | 0.34*** | −0.08 | 0.33*** | 0.41*** | 0.40*** | 0.27* | 0.72*** | 0.76*** | 0.65*** | 0.50*** | −0.31** | 0.95*** | 1.00 |
| Mean | 0.55 | 0.19 | 0.51 | 0.54 | 0.55 | 0.38 | 0.59 | 0.59 | 0.52 | 0.46 | 0.31 | 0.46 | 0.48 |
“mean” represents the average value of the absolute correlation coefficient of one index with others. *: corrected |$P<0.05$|, **: corrected |$P<0.01$|, ***: corrected |$P<0.001$|.
As shown in the first column of Fig. 2, CMY was highly correlated with all the functional indices across the cortex, especially with the static local functional features (ALFF, fALFF, and ReHo). Age, SA, and LGI were positively correlated with functional indices in a large portion of the cortex, except negatively correlated with d-ReHo around pre- and postcentral gyri, superior temporal cortex, precuneus, lateral and medial occipital cortex. CT was positively correlated with local functional features in the cingulate cortex, lateral posterior cortex, sensorimotor areas, primary visual cortex, and small parts of the lateral prefrontal cortex, and CT showed little correlation with connectivity-based features (FCS and stdFCS).
Correlation between age/structural and functional indices for each vertex. The significance was corrected using cluster-wise FDR correction (P < 0.05). The cold colors represent negative correlations while warm colors represent positive correlations. Blue and yellow represent that the absolute value of the correlation coefficient is larger than 0.5.
Regression analysis between structural and functional indices
Multiple regression models were built to explore the relationship between cortical structural and functional indices of neonates from the global, vertex, and region levels. For each level, the dependent variable was one of the functional indices, and the independent variables were made up of age, log(age) or age2, CT, SA, LGI, and CMY.
At the global level, CMY was positively linked to local functional features (maximal |$P$| = 0.0146, except d-ReHo,|$P$|= 0.8625, as shown in Table 2) and survive in backward elimination stepwise regression in models for fALFF, ReHo, d-ALFF and d-fALFF (maximal |$P$| = 0.0239, Supplementary Table S1). Models for fALFF and ReHo exhibited the best goodness-of-fit (adjusted R2 > 0.42) and the largest effect size (Cohen’s f2 > 0.7), while those for d-ReHo, FCS, and stdFCS exhibited the worst (adjusted R2 < 0.25 and Cohen’s f2 < 0.3). As covariates, BirthAge and weight were found significantly associated with ALFF, and BirthAge was significantly associated with FCS/stdFCS, so these covariates were controlled correspondingly in later analysis.
Regression results (standardized) of full models for the whole brain.
| Indices | AdjR2 | Effect size | CMY | CT | SA | LGI | Age | log(age) | Age2 |
|---|---|---|---|---|---|---|---|---|---|
| ALFF | 0.3006 | 0.4299 | 0.3550* | −0.0204 | −0.3525 | −0.2309 | 2.4048 | −1.6366 | – |
| fALFF | 0.4253 | 0.7401 | 0.4674*** | −0.0391 | −0.0425 | 0.0168 | 0.1011 | – | – |
| ReHo | 0.4574 | 0.8430 | 0.3705** | 0.0036 | −0.1597 | 0.2943 | 0.0166 | – | – |
| FCS | 0.1843 | 0.2259 | 0.0595 | −0.0621 | −0.0251 | 0.3167 | −0.1467 | – | – |
| d-ALFF | 0.3112 | 0.4519 | 0.4676** | 0.1229 | 0.1677 | −0.4843 | 4.4767** | – | −4.1449** |
| d-fALFF | 0.2938 | 0.4161 | 0.4005** | 0.1869* | 0.2940 | −0.6511* | 5.7292*** | – | −5.3126*** |
| d-ReHo | 0.1309 | 0.1507 | 0.0277 | 0.1838 | 0.1781 | −0.4856 | 3.5560* | – | −3.4904* |
| stdFCS | 0.2277 | 0.2949 | 0.0979 | −0.0134 | −0.1021 | 0.4850 | −0.2764 | – | – |
| Indices | AdjR2 | Effect size | CMY | CT | SA | LGI | Age | log(age) | Age2 |
|---|---|---|---|---|---|---|---|---|---|
| ALFF | 0.3006 | 0.4299 | 0.3550* | −0.0204 | −0.3525 | −0.2309 | 2.4048 | −1.6366 | – |
| fALFF | 0.4253 | 0.7401 | 0.4674*** | −0.0391 | −0.0425 | 0.0168 | 0.1011 | – | – |
| ReHo | 0.4574 | 0.8430 | 0.3705** | 0.0036 | −0.1597 | 0.2943 | 0.0166 | – | – |
| FCS | 0.1843 | 0.2259 | 0.0595 | −0.0621 | −0.0251 | 0.3167 | −0.1467 | – | – |
| d-ALFF | 0.3112 | 0.4519 | 0.4676** | 0.1229 | 0.1677 | −0.4843 | 4.4767** | – | −4.1449** |
| d-fALFF | 0.2938 | 0.4161 | 0.4005** | 0.1869* | 0.2940 | −0.6511* | 5.7292*** | – | −5.3126*** |
| d-ReHo | 0.1309 | 0.1507 | 0.0277 | 0.1838 | 0.1781 | −0.4856 | 3.5560* | – | −3.4904* |
| stdFCS | 0.2277 | 0.2949 | 0.0979 | −0.0134 | −0.1021 | 0.4850 | −0.2764 | – | – |
|$\mathrm{Adj}{R}^2$| : adjusted R square of the model; Effect size: Cohen’s f2(Cohen 2013). *: |$P<0.05$|, **: |$P<0.01$|, ***: |$P<0.001$|.
Regression results (standardized) of full models for the whole brain.
| Indices | AdjR2 | Effect size | CMY | CT | SA | LGI | Age | log(age) | Age2 |
|---|---|---|---|---|---|---|---|---|---|
| ALFF | 0.3006 | 0.4299 | 0.3550* | −0.0204 | −0.3525 | −0.2309 | 2.4048 | −1.6366 | – |
| fALFF | 0.4253 | 0.7401 | 0.4674*** | −0.0391 | −0.0425 | 0.0168 | 0.1011 | – | – |
| ReHo | 0.4574 | 0.8430 | 0.3705** | 0.0036 | −0.1597 | 0.2943 | 0.0166 | – | – |
| FCS | 0.1843 | 0.2259 | 0.0595 | −0.0621 | −0.0251 | 0.3167 | −0.1467 | – | – |
| d-ALFF | 0.3112 | 0.4519 | 0.4676** | 0.1229 | 0.1677 | −0.4843 | 4.4767** | – | −4.1449** |
| d-fALFF | 0.2938 | 0.4161 | 0.4005** | 0.1869* | 0.2940 | −0.6511* | 5.7292*** | – | −5.3126*** |
| d-ReHo | 0.1309 | 0.1507 | 0.0277 | 0.1838 | 0.1781 | −0.4856 | 3.5560* | – | −3.4904* |
| stdFCS | 0.2277 | 0.2949 | 0.0979 | −0.0134 | −0.1021 | 0.4850 | −0.2764 | – | – |
| Indices | AdjR2 | Effect size | CMY | CT | SA | LGI | Age | log(age) | Age2 |
|---|---|---|---|---|---|---|---|---|---|
| ALFF | 0.3006 | 0.4299 | 0.3550* | −0.0204 | −0.3525 | −0.2309 | 2.4048 | −1.6366 | – |
| fALFF | 0.4253 | 0.7401 | 0.4674*** | −0.0391 | −0.0425 | 0.0168 | 0.1011 | – | – |
| ReHo | 0.4574 | 0.8430 | 0.3705** | 0.0036 | −0.1597 | 0.2943 | 0.0166 | – | – |
| FCS | 0.1843 | 0.2259 | 0.0595 | −0.0621 | −0.0251 | 0.3167 | −0.1467 | – | – |
| d-ALFF | 0.3112 | 0.4519 | 0.4676** | 0.1229 | 0.1677 | −0.4843 | 4.4767** | – | −4.1449** |
| d-fALFF | 0.2938 | 0.4161 | 0.4005** | 0.1869* | 0.2940 | −0.6511* | 5.7292*** | – | −5.3126*** |
| d-ReHo | 0.1309 | 0.1507 | 0.0277 | 0.1838 | 0.1781 | −0.4856 | 3.5560* | – | −3.4904* |
| stdFCS | 0.2277 | 0.2949 | 0.0979 | −0.0134 | −0.1021 | 0.4850 | −0.2764 | – | – |
|$\mathrm{Adj}{R}^2$| : adjusted R square of the model; Effect size: Cohen’s f2(Cohen 2013). *: |$P<0.05$|, **: |$P<0.01$|, ***: |$P<0.001$|.
The vertex-wise regression results were shown in Fig. 3, CMY was positively associated with both static and dynamic local functional features across nearly the whole cortex. SA was negatively associated with ALFF in the left middle temporal gyrus. LGI was negatively associated with fALFF in certain sulci, such as the postcentral sulcus, superior and inferior temporal sulcus, and superior frontal sulcus, and was positively associated with fALFF and ReHo in the right precuneus and primary visual cortex. Age and log(age) were associated with ALFF in the bilateral cingulate cortex, the supramarginal gyrus, and the lateral prefrontal cortex. Results on the right hemisphere were presented in Supplementary Fig. S3.
Vertex-wise regression results between age, structural and functional indices (on the left hemisphere). The significance was corrected using cluster-wise FDR correction (P < 0.05). The cold colors represent -log(P) values corresponding to negative regression coefficients, while warm colors represent -log(P) values corresponding to positive regression coefficients. Blue and yellow represent the -log(P) value is larger than 10.
The brain was partitioned into 200 regions by applying the functional parcellation method that we proposed in a previous study (Zhao et al. 2020) to the neonatal brain, and the results of region-level regression analysis were largely similar to the vertex-level ones (details in Supplementary Fig. S4-5). CMY was positively associated with (d-)ALFF, (d-) fALFF, and ReHo across the whole cortex, and was associated with d-ReHo in smaller parts of the cortex. CT was positively correlated with dynamic local functional features in the circular sulcus of insula, small parts of the cingulate, temporal, and lateral prefrontal cortices, and the lateral and medial occipital cortices. For both vertex-wise and parcellation analysis, the best fitness (adjR2) was found in models for fALFF and ReHo (close to 0.3), while the worst was in models for FCS and d-ReHo (around 0.1, details in Supplementary Table S2).
Discussion
In the present study, we explored the relationship between brain structure and function development in neonates comprehensively via a variety of MRI indices. The results demonstrate a strong association between cortical myelination and local brain functional activity on multiple scales in both correlation and regression analyses. When we controlled the impact of common age dependence in models, cortical morphology is not as highly as CMY associated with local functional features in the large portion of the cortex.
Development of brain structure and function
The overall structural and functional indices of the neonatal brain (except d-ReHo) are increasing with age (Fig. 2 and Supplementary Fig. S2), suggesting the growing trend toward a bigger, more complex, more active, more regionally synchronized, and more temporally dynamic brain. Intrinsic resting-state brain activity is thought to predict subsequent task-evoked brain activations and behavioral performance (Fox et al. 2007; Zou et al. 2013), and higher local brain activity may increase the possibility of active brain responses to outer environmental stimuli, which is important for early development. Tight functional associations between one voxel and its neighbors guarantee effective collaboration to support specific functions, especially in primary sensory areas (Jiang et al. 2015), and local connectivity is rather essential before long-range connectivities become mature in the neonatal brain. Time-varying functional brain organization supports the need for rapid responses to the dynamically changing world, and higher temporal variability of the brain may reflect more flexibility of cognitive process in dealing with the environment, in line with a prior study that the overall functional dynamics were increasing in infancy (Wen et al. 2020).
Robust effects of CMY on local functional features
In addition to the main results, a series of validation analyses also demonstrated robust effects of CMY on local functional features of the neonatal brain. Global signal regression was applied to BOLD signals, and the regression results at both global and vertex levels were substantially similar to the main results (details in Supplementary Table S3–5 and Fig. S6–7). A smaller smooth kernel (4 mm) showed limited impact on the close association between CMY and local functional features at the global level (details in Supplementary Table S6). A larger (50 s) or smaller (30 s) window did not alter the relationship between CMY and dynamic ALFF/fALFF, while the models with a smaller window exhibited better fitting performance (details in Supplementary Table S7–8). Without preterm infants, similar to the main results, CMY was significantly associated with local functional features, but the model goodness-of-fit and effective size decreased, indicating the importance of including preterm infants and a broader age span (details in Supplementary Table S9). Since dHCP provided radiology scores for neonatal brain images, we included the scores as a covariate in the regression analysis, and the results still showed similar wide impact of CMY on local functional features (details in Supplementary Table S10). When subjects were randomly split into training and testing groups (details in Supplementary methods), the averaged correlation coefficient between the actual and predicted values of the testing datasets ranged from 0.28 to 0.65 within which the correlation coefficient for fALFF/ReHo was over 0.6 (details in Supplementary Table S11). Results of mediation analysis showed that CMY significantly mediated the relationship between age and fALFF/ReHo/d-ALFF (P < 0.01) (details in Supplementary methods and Fig. S8).
Relationship between CMY and local functional features
CMY is the most important and powerful index for functional brain features in this study. CMY was highly associated with local functional features on multiple scales. To note, different from the spatially varied coupling between structural and functional brain features discovered in adults (Huntenburg et al. 2017; Gu et al. 2021), the associations between CMY and brain functions were non-specific to areas or features of local functional activity, suggesting a fundamental role of CMY in early functional brain development.
The observed close association between CMY and functional brain activity is rooted in neurobiological bases. It is possible that increased conduction velocity of the nerve impulse along axons (Nave and Werner 2014) and reduced signal interference between axons (Hull et al. 2015) induced by better myelination makes signal transmission function more smoothly among microcircuits, which increases the probability of synchronized firing patterns of neuronal populations. Specifically, from the biological view (Almeida and Lyons 2017), the increased conduction velocity changes (i) the timing of arrival in postsynaptic neurons, (ii) the frequency of firing, (iii) the rate of depolarization and repolarization, which affects the probability of the generation of synchronous firing; and (iv) oscillations patterns among numerous neuronal populations, which further affect the macroscale functional activity and cognitive functions of the brain (Pajevic et al. 2014; Turner 2019). In addition, adding the spatial distribution of myelin to the modeling of functional connectivity can improve the goodness of fit (Demirtaş et al. 2019). Even without any predefined information, the cortical gradients of the modeling results of the large-scale circuit still appear strong relevance to the myelin content map (Wang et al. 2019). These studies above together with our results illustrate the close relationship between myelin and neuronal activity on both micro- and macroscale and the significant role of myelin in bridging between the two scales.
The influence of CT or SA on functional brain activity may be limited by the maturation of myelin in the cerebral cortex. Most neurons reach their destinations before birth and start to extend their axons and dendrites resulting in the explosion of synapse number (Keunen et al. 2017), which is reflected in increases in cortical thickness and expansions of surface area (Gilmore et al. 2018). Without myelin sheath around axons, the transmission of neural signals is inefficient when synapses grow dramatically in vertical (as revealed by CT) and horizontal directions (as revealed by SA). In other words, CMY may have a more direct impact on neural activity than CT or SA.
Relationship between CMY and connectivity-based features
Distinctions of measurement levels may lead to different results of structure–function correspondence. CMY predicts local connectivity (ReHo) much better than non-locals (FCS and stdFCS). Long-distance functional connections may be highly associated with white matter myelination, which are both very unmatured at birth (Dubois et al. 2014; Zhang et al. 2019), while local functional connection development might be easier and faster at this stage. Structural connectivity architecture was found capable to predict static functional connectivity (Honey et al. 2009) or shape dynamic states of high modularity (Liégeois et al. 2016). Pair-wise or state-wise exploration of the relationship between CMY and functional connectivity might provide more insights.
Relationship between cortical morphology and functional features
In comparison with CMY, cortical morphology is not that closely associated with local brain functional activity in neonates, which is inconsistent with the results in adults (Jiang et al. 2015; Qing and Gong 2016). A previous study modeled the relationship between ReHo and several morphological features such as CT, SA, and LGI, and found large-scale structure–function correlations in adults (Jiang et al. 2015). Our goodness-of-fit (0.25–0.28 for the average adjR2) of models for ReHo outperforms the previous study (<0.12) (Jiang et al. 2015), indicating the advantage of including CMY in regression models for the neonatal brain. The correlation between morphological features and between local functional features is stronger than the cortical morphology-functions correlation, consistent with (Yang et al. 2016). Cortical morphological features are widely used to explore their association with the development of cognitive functions and brain activity (Shaw et al. 2006; Geng et al. 2017; Girault et al. 2020). Our results indicate the essential role of CMY in functional development in the local scales and suggest more attention to the CMY index when investigating the biological basis of functional brain activity, cognitive ability, behavioral performance, and neurodevelopmental disorders.
Limitations
Several limitations should be noted in this study. First, the number of very young (|$\le$| 36 weeks) or old (44 weeks) newborns is disproportionately smaller than other ages, which may affect the accuracy of the model results. Second, the functional periods that we used were shorter than those in other dynamic functional connectivity studies, though the number of volumes seemed sufficient (Allen et al. 2014; Yin et al. 2020). The short period of functional data is due to the special case of infant brain scanning, in which big head movement occurs occasionally during infants’ natural sleep. Third, the frequency band (0.01–0.1 Hz) was used in the computation of ALFF and fALFF, which may not be ideally suitable for neonates (Alcauter et al. 2015), and results using different frequencies would be better compared. Fourth, infants were scanned during sleep, and better structure–function coupling was found in sleep or anesthesia than being awake (Barttfeld et al. 2015; Tagliazucchi et al. 2016). The conscious state of sleep might be a confounder for the results, but we could still claim a more dominant impact of cortical myelination on functional activity than other variables. To note, the effect of age on functional activity might be limited by the narrow age range (29–44 weeks) in this study. Finally, multiple regression models were constructed in this study, and it is acknowledged that there may be more elaborate models to better delineate the relationship between neonatal brain structure and function. To figure out whether the strong association between CMY and brain activity is specialized for neonates or not, studies in adults and children are expected to conduct.
Conclusion
The present study demonstrates the close relationship between cortical myelination and local functional features of the neonatal brain. The results highlight the fundamentally important role of cortical myelination compared to age or cortical morphology in early functional brain development. Our study provides valuable insights into understanding the brain structure–function relationship at birth, suggests more attention to myelination in studies of early brain development, and helps bridge between anatomical abnormalities and neuropsychiatry disorders.
Funding
This work was supported by the Key Realm R&D Program of Guangdong Province (2019B030335001), NFSC (National Natural Science Foundation of China) (Grant No. 61403148). Neonatal data were provided by the developing Human Connectome Project, KCL-Imperial- Oxford Consortium funded by the European Research Council under the European Union Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement no. [319456].
Conflict of interest statement: None declared.
