Examining the relationships between cortical maturation and white matter myelination throughout early childhood

Cortical development and white matter myelination are hallmark processes of infant and child neurodevelopment, and play a central role in the evolution of cognitive and behavioral functioning. Non-invasive magnetic resonance imaging (MRI) has been used to independently track these microstructural and morphological changes in vivo, however few studies have investigated the relationship between them despite their concurrency in the developing brain. Further, because measures of cortical morphology rely on underlying gray–white matter tissue contrast, which itself is a function of white matter myelination, it is unclear if contrast-based measures of cortical development accurately reflect cortical architecture, or if they merely represent adjacent white matter maturation. This may be particularly true in young children, in whom brain structure is rapidly maturing. Here for the first time, we investigate the dynamic relationship between cortical and white matter development across early childhood, from 1 to 6 years. We present measurements of cortical thickness with respect to cortical and adjacent myelin water fraction (MWF) in 33 bilateral cortical regions. Significant results in only 14 of 66 (21%) cortical regions suggest that cortical thickness measures are not heavily driven by changes in adjacent white matter, and that brain imaging studies of cortical and white matter maturation reflect distinct, but complimentary, neurodevelopmental processes.


Introduction
Two important neurodevelopmental processes that occur throughout infancy and early childhood are the maturation of myelinated white matter and the development of the cerebral cortex. The formation of the lipid bilayer myelin sheath around neuronal axons (myelination) is essential for the rapid brain messaging required for higher order behavioral and cognitive functioning. Brain disconnectivity resulting from aberrant or insufficient development of the myelin sheath may underlie a number of neuropsychiatric disorders, including autism and attention deficit hyperactivity disorder (Belmonte et al., 2004;Krain and Castellanos, 2006;Konrad and Eickhoff, 2010;Xiao et al., 2014). Measures of cortical development, including changes in thickness, surface area, gyrification, volume, and gray matter myelination, have also been investigated in the context of developmental and psychiatric disorders (Courchesne et al., 2007;Hazlett et al., 2012), and in relation to cognitive performance (Shaw et al., 2007.
Advances in magnetic resonance imaging (MRI) have allowed for the in vivo investigation of myelination and cortical maturation both across development and in association with cognitive and behavioral development. Multicomponent relaxometry (MCR) techniques, such as mcDESPOT (multicomponent driven equilibrium single pulse observation of T 1 and T 2 ) (Deoni et al., 2008), enable the visualization and quantification of a surrogate measure of myelin content, termed the myelin water fraction (MWF). MCR decomposes the measured MRI signal into the contributions of signal signatures associated with differing microanatomical water compartments. In the brain, three distinct water pools are commonly observed, corresponding to the free intra and extra-axonal water, the CSF water, and the water trapped between lipid bilayers of the myelin sheath (MacKay et al., 2006). Quantification of the myelin-associated signal, the MWF, is a useful metric for tracking white matter maturation (Deoni et al., 2012;Dean et al., 2015) and its relationship to cognitive development (O'Muircheartaigh et al., 2013 in the developing brain. Most commonly measured through Freesurfer segmentation (Fischl, 2012) of the cortical ribbon from a T 1 -weighted MR image, cortical thickness is an oft used metric for tracking synaptic density and cortical maturation. While Freesurfer analysis is not formally recommended for use in children under 4 years of age, it has been utilized in investigations of cortical development in infants and toddlers as young as 12 months (Lowe et al., 2012;Travis et al., 2014). Accurate and reproducible delineation of cortical gray matter from underlying and adjacent white matter is a prerequisite for calculating cortical thickness. In ours (Deoni et al., 2015), and others (Lyall et al., 2014) experience, inaccuracies in cortical segmentation can be attributed to insufficient gray-white matter contrast in children under 12 months.
While myelination and cortical development do not occur independently, with both processes occurring symbiotically during the development of neural systems, few studies have sought to investigate the relationship between them. Further, since accurate cortical thickness measures necessitate strong gray-white matter image contrast, which itself is a function of white matter myelination, it is unclear if measures of cortical thickness in early childhood reflect cortical architecture or adjacent white matter maturation. In this work, we aim to directly examine the relationships between cortical thickness and white matter myelination in a large cohort of 134 typically-developing children between 1 and 6 years of age. We measured cortical thickness and calculated the MWF within directly adjacent white matter in 33 bilateral cortical regions. Our results show that cortical thickness changes are not fully explained by MWF changes alone, suggesting that Freesurfer cortical thickness values and MWF are measuring distinct and complementary processes of neurodevelopment.

Study design and participants
Data from 134 (58 female) healthy and typically-developing children approximately 1 to 6 years of age (363 to 2198 days corrected to a 40-week gestation) were used in this study. These children were recruited as part of an ongoing longitudinal investigation of white matter maturation in relation to behavioral development in infancy and early childhood (Deoni et al., 2012). Full demographic information is provided in Table 1. A total of 177 scans were performed, with 36 children scanned at least twice and 7 children scanned three times. The average time between repeat scans was approximately one year ( Fig. 1). Inclusion criteria consisted of: birth between 37 and 42 weeks gestation; no abnormalities present on fetal ultrasound; no delivery complications (i.e. no visits to the neonatal intensive care unit); APGAR score of 8 or higher; no in utero exposure to illicit drugs or alcohol; no pregnancy complications (i.e. preeclampsia); no familial history of learning disability, behavioral or psychiatric disorder; and no reported neurological events or disorders in the infant such as head trauma or epilepsy. Child, sibling, and parent medical histories were collected as a supplement to parental interviews conducted at the time of study enrollment. Written informed consent was obtained from the parent(s) or legal guardian of each participating child, and all experimentation was performed under the approval of the Brown University Institutional Review Board.

Image acquisition
To measure MWF, whole-brain mcDESPOT data were acquired using age-optimized imaging protocols described previously (Deoni et al., 2012) and summarized in Table 2. All imaging was performed on a 3T Siemens Tim Trio scanner with a 12-channel head RF coil array. The data set for each child includes 8 T 1 -weighted spoiled gradient echo (SPGR) images, 2 inversion-prepared SPGR images (IR-SPGR), and 2 sets of T 1 /T 2 -weighted steady-state free precession (bSSFP) images, each acquired with a differing radio-frequency phase-cycling pattern (Deoni, 2011). High resolution volumetric T 1 -weighted MP-RAGE data were also acquired for cortical morphometry analysis.
Children under the age of four were imaged during natural (non-sedated) sleep, while children over four were imaged while watching a favorite movie or TV show . To attenuate noise levels in the scanner, and keep the youngest participants asleep for the duration of the session, peak gradient amplitudes and slew rates were reduced to 25 mT/m/s, foam inserts (Quiet Barrier HD Composite, UltraBarrier USA) were secured inside the scanner bore, and headphones (MR Confon, Germany) were positioned to cover the ears. To limit the possibility of movement during the scan, all children were swaddled in ageappropriate MedVac vacuum immobilization bags (CFI Medical Solutions, USA) and their heads were kept in place with foam pads. A research assistant watched over infants from inside the scanner room, and additional monitoring was possible using a pediatric pulseoximetry system and infrared camera. During acquisition, image data was evaluated for motion artifacts including blurring and ghosting. Presentation of these artifacts on an image necessitated repeated acquisition of that image using the original FOV positioning and sequence parameters until higher quality data was obtained. These motion-free images were then incorporated into the child's data set as replacements for artifact-laden images prior to image processing .

Image analysis
Generation of myelin water fraction maps-Following image acquisition, the mcDESPOT data from each child were linearly co-registered to account for subtle intra-scan motion (Jenkinson et al., 2002), non-brain signal was removed (Smith, 2002), B 0 and B 1 field calibration maps were calculated (Deoni, 2011), and voxel-wise T 1 and MWF maps were calculated using an iterative stochastic region contraction approach (Deoni and Kolind, 2014).
Delineation of cortical regions-Low frequency intensity inhomogeneities were removed from the high resolution T 1 -weighted images using Advanced Normalization Tools (ANTs) nonparametric non-uniform normalization (N3) bias correction. Freesurfer (Fischl, 2012) was then used to demarcate the cortex into 33 distinct regions per hemisphere and segment the cortical ribbon for cortical thickness calculations. At each stage in the Freesurfer processing pipeline, the MP-RAGE T 1 -weighted images were visually inspected and, if needed, manually edited and corrected. This included inspecting data for poor skullstripping, the additional use of gcut (http://freesurfer.net/fswiki/FsTutorial/ SkullStripFix_freeview) and, in extreme cases, manual removal of remaining dura, eye, and other non-brain signal. Mean cortical thickness values were obtained for each region.
Calculation of adjacent white matter MWF-To obtain measures of the mean adjacent white matter MWF for each cortical region, each region mask was first blurred with a 2 dimensional Gaussian kernel with a 4 mm FWHM. Non-white matter signal was removed by subtracting the cortex and any other non-brain regions. This blurred mask was then superimposed on the aligned MWF map and the mean MWF value was calculated.
Calculation of cortical MWF-Mean MWF values were obtained by co-registering each child's MWF map to their high resolution T 1 weighted image and superimposing each delineated region onto this registered map. Non-linear registration (Avants et al., 2011) was performed first using the high flip angle T 1 weighted SPGR image acquired as part of the mcDESPOT protocol, with the transformation matrix subsequently applied to the MWF map. Mean and standard deviation MWF values were calculated for each region, in each hemisphere, for each child.
Similar analysis was performed for the quantitative T 1 maps calculated as part of the mcDESPOT processing, with mean T 1 values obtained for each of the 66 cortical and adjacent white matter regions.
A visual overview of each image analysis step is provided in Fig. 2.

Region development trajectories
For each region, plots of (1) cortical thickness and adjacent white matter MWF vs. age; (2) cortical thickness and cortical MWF vs. age; and (3) cortical and adjacent white matter MWF vs. age were generated. We fit continuous logarithmic functions to the MWF and cortical thickness data in these plots, verifying the most appropriate fit to the cortical thickness data using the Bayesian Information Criterion (BIC).
To examine the relationships between measures, we first removed the effect of age on each measure by fitting the appropriate trend and subtracting it from the raw values (i.e., calculating the residuals). We then calculated the Pearson product-moment correlation coefficient between these residuals for (1) cortical thickness vs. adjacent white matter MWF; (2) cortical thickness vs. cortical MWF; and (3) cortical vs. adjacent white matter MWF. The first of these analyses examined whether the maturation of adjacent myelin is a significant driver for cortical growth; the second sought to determine if changes in cortical myelin content are a significant driver for cortical growth; and the third explored the strength of the association between cortical and adjacent white matter development. We accounted for multiple comparisons using the Holm-Bonferroni method with an alpha of 0.05 (corrected for 33 bilateral region comparisons).
Using the average longitudinal relaxation (T 1 ) times for each of the 66 cortical and adjacent white matter regions, we also calculated the expected ideal (i.e., without contaminating proton density effects) T 1 -weighted signal as (1-2e -TI/T1 ), with TI = 950 ms, matching the TI of the acquired MP-RAGE data. Gray-adjacent white matter T 1 contrast was then calculated for each of the 66 regions, and this contrast plotted against the region's mean cortical thickness. This analysis was performed to determine the sensitivity of, and relationship between, cortical thickness measures and ideal image contrast. As before, the Pearson's r was calculated and assessed for significance using the Holm-Bonferroni method with an alpha of 0.05. Fig. 3 shows raw data plots and superimposed growth models for (1) cortical thickness and adjacent white matter MWF vs. age; (2) cortical thickness and cortical MWF vs. age; and (3) cortical and adjacent white matter MWF vs. age for a representative subset of slow, moderate, and fast developing cortical regions. In agreement with prior data obtained by our group (Deoni et al., 2012(Deoni et al., , 2015, we model the development of both cortical and adjacent white matter MWF using an increasing logarithmic function. The fit curve equations in Table 3 reveal a range in the logarithmic slope (i.e. the rate of MWF development) across cortical regions. Absolute cortical myelin content values are, overall, lower compared to adjacent white matter, as expected.

Results
To investigate how cortical thickness changes with age, we fit logarithmic, quadratic, and linear growth models to the data and compared them using the BIC. While some regions were more appropriately characterized by linear or quadratic fits, the majority of regions follow a logarithmic trajectory (Table 4). To us, these results justified global logarithmic modeling of cortical thickness for all subsequent analyses. Apart from 5 regions (bilateral entorhinal, right parahippocampal, and bilateral temporal pole), measures of cortical thickness follow a decreasing trajectory with age.
The relationships between cortical thickness, adjacent white matter MWF, and cortical MWF are shown in Fig. 4, which contains plots of the residuals for each measure against the others, and Table 5, which details quantitative results from correlation analyses. By examining the residuals (calculated by subtracting the logarithmic model predictions from measured values), we removed the effect of age from the data. Comparing changes in cortical thickness with those of adjacent white matter reveals a significant (p < 0.05 corrected for multiple comparisons) negative relationship between these processes in 10 of 66 regions (Pearson's r range: −0.374 to −0.252), including the inferior parietal, supramarginal, rostral middle frontal, and superior frontal regions bilaterally. A significant positive relationship was found in 4 of 66 regions (Pearson's r range: 0.249 to 0.29), including the right cuneus, right lingual, and bilateral transverse temporal regions.
In 16 of 66 total regions, we found a statistically significant (p < 0.05 corrected for multiple comparisons) negative relationship between cortical thickness and cortical MWF (i.e. greater thickness is associated with lower MWF). Correlation coefficients for these relationships range from −0.248 to −0.474. Bilateral significance in the pars triangularis, caudal middle frontal, middle temporal, inferior parietal, inferior temporal, and supramarginal regions accounts for 12 of these results. The remaining 4 significant relationships are found in the right postcentral, left rostral middle frontal, left superior temporal, and left parsopercularis regions. Finally, cortical MWF and white matter MWF show significant positive relationships (Pearson's r range: 0.209 to 0.742) in 63 of 66 regions. Here, the only non-significant regions are the left entorhinal cortex and bilateral temporal pole.
Significant relationships (p < 0.05 corrected for multiple comparisons) between T 1 contrast and cortical thickness exist in 10 of 66 regions, including both hemispheres of the inferior parietal, middle temporal, and pars orbitalis regions. Globally, T 1 contrast varies only subtly between 1 and 6 years of age. Similarly to Figs. 3 and 4, Fig. 5 highlights this analysis for the left hemisphere superior parietal, supramarginal, and middle temporal regions.

Discussion
In this work, we have investigated the dynamic relationship between cortical development and white matter maturation using quantitative high resolution and MWF imaging for the first time. In a large cohort of 134 (58 female) healthy and typically-developing children, we show that cortical thickness, cortical myelin, and adjacent white matter myelin each follow logarithmic development trajectories. Myelin trajectories presented here are consistent with our prior investigations (Deoni et al., 2012(Deoni et al., , 2015 and cortical thickness trajectories were chosen through BIC analysis. Visual inspection of residual plots revealed approximately normal distributions of points around the origin, providing additional evidence in favor of these models. In 61 of the 66 regions examined, cortical thickness is found to decrease logarithmically from 1-6 years of age. Prior studies have demonstrated early expansions in cortical development from birth to 1 year of age, and region-specific cortical thinning from 1-2 years of age (Lyall et al., 2014). Our results reveal that visual, motor, and somatosensory areas appear to have faster rates of cortical thinning compared to frontal and association regions, although future studies are needed to verify these trends. These changes in cortical thickness occur simultaneously with logarithmic increases in both cortical MWF and adjacent white matter MWF. Prior work has revealed a similar relationship between cortical thinning and brain growth in children between the ages of 5 and 11 (Sowell et al., 2004). While the authors suggest that this is perhaps due to increased cortical myelination in lower cortical layers, our analysis presents a more complex picture of cortical development and myelination during the first few years of life.
In particular, we find that cortical development is significantly correlated with both cortical white matter and adjacent white matter maturation in relatively few regions. These relationships are not concentrated in one brain area but instead can be found across the brain in regions that differ in rates of cortical thinning. Notably, these regions include latermyelinating frontal and association regions such as the inferior parietal, supramarginal, rostral middle frontal, and caudal middle frontal regions. When considering all brain regions, however, the square of the correlation coefficient, r 2 , between cortical thickness and adjacent white matter MWF does not exceed 0.14. In this case, the linear regression model used to illustrate the relationship between the residual measurements only accounts for 14% of the variability in the data, suggesting that cortical thickness and adjacent white matter MWF are not merely proxies for one another. Similarly, at most only 23% of the data is accounted for when explaining cortical thickness changes with respect to cortical MWF maturation. Combining these results with our observation that the majority of cortical regions show non-significant correlations suggests that measures of cortical thickness and MWF are complimentary, but do not characterize identical underlying processes.
Further support for this claim comes from an analysis of the relationship between cortical thickness and gray-white matter T 1 contrast. Across early childhood (specifically up to 5 years of age), white matter myelination advances in a caudal-cranial, posterior-anterior pattern. In combination with changes in fiber density and coherence, compartmentalization of free water, and changes in macromolecule, protein, lipid, and cholesterol content, this maturation results in significant reductions in white matter relaxation parameters (T 1 and T 2 ). Analogously, changes in synaptic density and cortical architecture also result in widespread reductions in cortical T 1 and T 2 (Deoni et al., 2015). Combined, these changing MRI parameters yield a maturing gray-white matter tissue contrast that gradually takes on an adult-level appearance (Barkovich et al., 1988;Paus et al., 2001). Knowing that myelination contributes to changes in T 1 , evidence of a relationship between T 1 contrast and cortical thickness could be suggestive of a developmental connection between myelination and cortical thickness. However, our analysis shows little evidence of such a relationship. This lends support to measures of cortical thickness being independent of white matter MWF. Looking further at T 1 contrast values over time, we also see that adult levels of contrast are established and relatively stable by 1 year of age, compared to the increasing logarithmic trajectory of both cortical and white matter MWF from 1 to 6 years of age. This may suggest that myelination is not the primary driver of T 1 contrast within this age range, a conclusion supported by prior null findings of a MWF-T 1 relationship in white matter across childhood (Deoni et al., 2012;Harkins et al., 2015).
A potential methodological concern with this work lies in the relatively low resolution of the Freesurfer (Fischl, 2012) segmented cortical regions. To ensure accurate parcellation, images were visually inspected at each stage in the processing pipeline. Children under the age of 1 were also excluded from this work due to insufficient gray-white matter contrast observed in this age range. While lower than the recommended 1 mm 3 isotropic resolution for adult studies, the (1.2 × 1.2 × 1.2) mm 3 spatial resolution of our T 1 w images either meets or exceeds resolutions used in prior pediatric neuroimaging studies (Deoni et al., 2015;Shaw et al., 2012).
Multicomponent relaxometry techniques, such as mcDESPOT, are specific to early myelin development (Deoni et al., 2008(Deoni et al., , 2012(Deoni et al., , 2013. White matter microstructural changes, however, extend beyond myelination and encompass changes in axon number and density. Diffusion tensor imaging (DTI) can provide insight into these additional neuroanatomical measures, but sacrifices myelination specificity (Mädler et al., 2008;Yoshida et al., 2013). Future studies using diffusion tensor imaging (DTI) in combination with mcDESPOT are needed to gain a more comprehensive understanding of early white matter development in this age range.
While this work highlights primarily non-significant relationships between measures of cortical development and white matter maturation, a temporal offset may exist between these processes that was not considered here. Prior work has shown that over time, trajectories of cortical thickness changes are regionally-dependent and are associated with cognitive development and outcome (Shaw et al., 2006). Further analysis of both morphological and behavioral measures is necessary to examine whether early changes in cortical thickness may predict later changes in MWF, or vice versa. Gender is another factor that was not considered here. We have previously shown (Deoni et al., 2015) no significant evidence for sexual dimorphism in cortical MWF and T 1 development trajectories or mean values from 1 to 6 years of age. While sex-specific differences in the magnitude of cortical thickness have been observed from age 6 into adulthood, rate of cortical thickness change does not show gender influences in this period (Raznahan et al., 2011). Gender differences in cortical development and white matter maturation relationships may, therefore, be best investigated in late childhood and early adolescence, which is beyond the scope of this work.

Conclusions
Our results show that changes in cortical thickness from 1-6 years of age are non-linear and largely independent of both cortical and adjacent white matter maturation. These findings raise questions about the degree to which other cortical measures explain the relationship between cortical and white matter development. While further investigation is needed to determine if the regional variation in cortical thickness shown here can be linked to cognitive and behavioral outcomes, our results fill in the knowledge gap on cortical and white matter development trajectories and their relationship to one another in early childhood. Age distribution (corrected to a 40-week gestation) of study cohort with females in green and males in blue. Individual scans are denoted by an asterisk, with dashed lines connecting repeated measurements from the same child.

Fig. 2.
Image analysis steps. (a) MWF maps were calculated using an iterative stochastic region contraction approach. Freesurfer (Fischl, 2012) was used to demarcate the cortex (b) into 33 distinct regions per hemisphere and segment the cortical ribbon (c) for cortical thickness calculations. Freesurfer-derived cortical regions (d) were blurred with a 4 mm FWHM Gaussian kernel (e), and then gray matter and non-brain portions were removed (f). The final mask was then superimposed on to the co-registered MWF image (g) and mean white matter MWF was calculated.

Fig. 3.
Maturation profile comparisons between cortical thickness, cortical MWF, and adjacent white matter MWF in the left hemisphere of three cortical regions that vary in rate of cortical thinning. Similar trends are observed in the right hemisphere of these regions and in both hemispheres of the remaining 30 bilateral regions not pictured.

Fig. 4.
Comparisons between cortical thickness, cortical MWF, and adjacent white matter MWF residual values obtained by subtracting the logarithmic model predictions from measured values shown in Fig. 2. Asterisks denote a statistically significant (p < 0.05 corrected for multiple comparisons) relationship between the two measurements shown in a given plot.