The Concurrence of Cortical Surface Area Expansion and White Matter Myelination in Human Brain Development.

The human brain undergoes dramatic structural changes during childhood that co-occur with behavioral development. These age-related changes are documented for the brain's gray matter and white matter. However, their interrelation is largely unknown. In this study, we investigated age-related effects in cortical thickness (CT) and in cortical surface area (SA) as parts of the gray matter volume as well as age effects in T1 relaxation times in the white matter. Data from N = 170 children between the ages of 3 and 7 years contributed to the sample. We found a high spatial overlap of age-related correlations between SA and T1 relaxation times of the corresponding white matter connections, but no such relation between SA and CT. These results indicate that during childhood the developmental expansion of the cortical surface goes hand-in-hand with age-related increase of white matter fiber connections terminating in the cortical surface.

Introduction

Structural brain changes in human development are extensive and multifaceted. They include changes in cortical gray matter and also in the underlying white matter. These brain structural changes have been shown to correlate with the development of cognitive functions. For instance, such effects have been reported for the local increase of the white matter organization (Yeatman et al. 2012; Skeide and Friederici 2016; Grosse Wiesmann et al. 2017) and the local decrease of gray matter volume (GMV; Richardson et al. 2009).

The cortical gray matter contains most of the neuronal brain functions. These neurons are interconnected with each other and form short-range as well as also long-range network connections that are formed by the nerve fibers of the white matter. In early childhood, a marked increase in GMV is apparent, followed by a decrease after the age of 7 years (Giedd et al. 1999; Sowell et al. 2003).

Technically, GMV is defined as the product of cortical surface area (SA), i.e., the area covered by the cerebral cortex, on the one hand, and cortical thickness (CT), i.e., the combined thickness of the six layers of the human neocortex on the other. GMV appears as a rather gross measure, for it cannot distinguish concurrent changes in SA and CT. For this reason, a deeper insight into the gray matter composition may be gained by shifting the focus from volume measures to measures of CT and SA directly. SA, like CT, is a highly heritable trait (Eyler et al. 2011). Importantly, SA is independent of CT at the regional level as well as at the whole-brain level, and their genetic determinants are largely non-overlapping (Panizzon et al. 2009; Winkler et al. 2010). In fact, the two measures show little co-variation (Pakkenberg and Gundersen 1997; Im et al. 2008). As a consequence, each of these two measures may provide different pieces of information about developmental processes.

It has been shown that SA and CT follow distinct developmental trajectories (Hutton et al. 2009; Li et al. 2013; Wierenga et al. 2014; Lyall et al. 2015). While CT changes at a rather fast pace in the first two years of life, reaching 97% of adult CT by the age of 2 years (Lyall et al. 2015), SA shows much slower changes. SA increases to 69% of adult size at the age of 2 years (Lyall et al. 2015), and reaches adult levels only at about 8 years of age (Raznahan et al. 2011). The observed different trajectories of cortical change suggest that the CT and SA measure different aspects of the cortical development, respectively. Recent data have shown that high-expanding cortices are implicated in higher cognitive functions that show extensive development during childhood (Amlien et al. 2014; Fjell et al. 2015). A description of the mechanisms responsible for SA expansion during development, however, is not determined yet.

In his radial unit hypothesis, Rakic postulated that the size of cortical SA is driven by the number of cells within columnar units, hence neurogenesis is accounting for area growth (Rakic 1988). However, processes of neurogenesis and migration in the neocortex are mainly completed within the first weeks of postnatal life (Rakic 2009; Clowry et al. 2010; Hill et al. 2010). Thus, additional mechanisms are needed to explain SA increase as observed in further brain development. Cellular processes such as synaptogenesis, gliogenesis, dendritic arborization, and intra-cortical myelination are likely candidates that help shape cortical SA (Mrzljak et al. 1990; Rakic et al. 1994; Hill et al. 2010; Petanjek et al. 2011).

The developmental increase in GMV early in life is complemented by a volume increase in white matter that displays a steady correlation with age across childhood, adolescence, and early adulthood (Paus 2010). Several candidate mechanisms have been assumed to contribute to these age effects. Myelin sheath thickness is known to grow with age (Benes et al. 1994) and to continue until early adulthood (Yakovlev and Lecours 1967). Increase of white matter volume might also be dependent on changes in axonal caliber, which increases postnatally late childhood (Schröder et al. 1988). However, the relative contribution of myelin sheath and axonal caliber to total white matter volume might be difficult to disentangle since myelin-associated proteins modulate the caliber of myelinated axons (Yin et al. 1998).

Although both cortical SA and white matter volume show a more prolonged developmental trajectory compared with the rapid increase in CT, the relationship between cortical SA and the properties of the connected white matter axons has so far only marginally been considered (Seldon 2005). Given the similarities in the developmental trajectories of cortical SA changes and white matter maturation, a possible relation between these two developmental patterns can be assumed. In the present study we investigate whether cortical morphometric properties and myelin content of white matter, reflected by qT1, share a link. Here, we hypothesize that cortical SA expansion during development is related to the myelination of the axons originating from or terminating in the cortex. In order to test this, we analyzed gray matter morphometry and white matter composition and their relation to age in young childhood at 3–7 years. We chose this age group for two reasons: because gray and white matter of the brain has been shown to develop markedly at this age (Gogtay et al. 2004; Giedd and Rapoport 2010), and because their relation so far has not been considered for children younger than 5 years of age. To this end, we obtained cortical SA and thickness data in a group of N = 170 children to describe age effects in the morphology of gray matter, and T1 relaxation time data to measure age effects in white matter myelin content. We expected, first of all, to support previous findings showing SA expansion and CT decrease in brain maturation (Amlien et al. 2014). Importantly, we hypothesized that regional cortical surface expansion during childhood is linked to the maturation of white matter fibers that originate or terminate from the very region. If this was the case, white matter connections tracked from the gray/white matter interface, which show substantial myelin increase, should reveal cortical regions with prominent SA expansion. Such an analysis should result in a consistent pattern of surface expansion and myelination of connected pathways. To the best of our knowledge, this is the first study that puts in relation cerebral morphometric age-related effects with the myelin content of the deep white matter during early childhood.

Methods

Participants

MRI data were selected from a larger database collected at the Institute from several studies. There were no differences in scanner software, procedures, or experimenters involved in these studies. The database contained 333 sets of MRI data from 215 participants, acquired within 2 years. From this database, scans of children participating in training studies were excluded (48 datasets). Datasets displaying motion artifacts in either of the scans needed for the analysis were excluded (91 datasets). Twelve datasets were excluded for diagnostic reasons. Datasets belonging to the same participants (12 datasets) were excluded in order to balance age distribution. This resulted in a set of MRI data from 170 participants (88 girls, aged from 3.0 to 7.0 years, mean age 5.05 ± 1.15). This age range allowed observing marked age correlations in SA and T1 relaxation times, while still being restricted enough to approximate the global developmental patterns of these measures to a linear model. About one week before scanning, participants were acclimated to the MRI environment through a phantom scanner available at the Institute. After a short play session, participants observed a mock MRI scan being performed on a large plush doll, familiarizing with the equipment (earplugs, headphones, 3D goggles), the noise, and the procedure. After familiarization, participants underwent a mock scan themselves. During both the mock and the actual scan later, participants were allowed to watch an animated movie of their choice through the use of headphones and 3D goggles. All equipment was identical between mock and actual scans. Written informed parental consent was obtained for all participants prior to the experiment in accordance with the ethical approval from the University of Leipzig’s Ethics Committee.

Data Acquisition

High-resolution 3D T1-weighted images and diffusion-weighted data were acquired using a Siemens 3T TIM Trio scanner with a 12 channel array head coil. A high-resolution 3D T1-weighted image (MPRAGE) was acquired (inversion time, TI = 740 ms; repetition time, TR = 1480 ms; echo time TE = 3.46 ms; excitation flip angle α = 10°; image matrix = 256 × 240; field-of-view, FOV = 256 mm × 240 mm; slab thickness = 192 mm; 128 slices; sagittal orientation; spatial resolution = 1 mm × 1 mm × 1.5 mm; total acquisition time, TA = 5:57 min; no GRAPPA; bandwidth 190 Hz/Px).

Furthermore, in order to obtain reliable T1 maps, a 3D MP2RAGE sequence was acquired (Marques et al. 2010) (TI1 = 700 ms; α1 = 4°; TI2 = 2500 ms; α2 = 5°; repetition time of the total sequence cycle, TR = 5000 ms; TE = 2.82 ms; image matrix = 192 × 168; FOV = 250 mm × 219 mm; slab thickness = 188 mm; 144 slices; sagittal orientation; spatial resolution = 1.3 mm × 1.3 mm × 1.3 mm; TA = 6:22 min; generalized auto-calibrating partially parallel acquisitions, GRAPPA 3; bandwidth 606 Hz/Px). Diffusion-weighted images were acquired with twice-refocused spin echo echo-planar-imaging sequence (Parameters: TE = 83 ms; TR = 8 s; 100 × 100 image matrix; FOV = 186 × 186 mm2; 66 axial slices -no gap; resolution: 1.86 × 1.86 × 1.9 mm3; GRAPPA 2; fat saturation; partial Fourier factor 6/8; TA 9:20 min; bandwidth 1428 Hz/Px; echo spacing 0.78 ms). Two sets of diffusion-weighted images were acquired. In the first, acquisition was acquired along an anterior-to-posterior phase-encoding direction. Diffusion weighting was isotropically distributed along 60 diffusion-encoding gradient directions with a b-value of 1000 s/mm2. Additionally, seven images with no diffusion weighting (b0) were acquired initially and interleaved after each block of 10 diffusion-weighted images as anatomical reference for off-line motion correction. The second dataset was acquired along the posterior-to-anterior phase-encoding direction, and contained one b0 image along with one b1000 diffusion-weighted volume (TA 0:42). The second dataset was used to correct for the geometrical distortions caused by inhomogeneity of the static magnetic field.

Template Creation

Individual T1-weighted images were aligned to the MNI template using a rigid-body registration. Bias correction was performed using the N4 algorithm (Tustison et al. 2010). A common template has been created from the T1-weighted datasets using ANTS’s (http://stnava.github.io/ANTs/) buildtemplateparallel.sh script. The four default nonlinear registration iterations were used, plus an initial rigid-body registration. In order to reduce high frequency deformations, a symmetric diffeomorphic transformation was used with the following parameters: total variance in voxel space = 3; gradient step = 0.25. Registration parameters were tested on a sample of ten pediatric T1-weighted images registered to a fixed image. The parameters used in template creation were selected based on a combined assessment of visual inspection and cross-correlation values between the fixed and the moving images

The resulting group template was processed in Freesurfer to reconstruct the cortical surfaces. Surfaces were processed following the standard Freesurfer’s pipeline. Segmentations were inspected manually for accuracy, and corrected when needed. Errors usually pertained to misclassification of few voxels near membranes or adjacent vessels. Cortical and pial surfaces were inspected for errors, and manual correction through control points was performed. The automated Freesurfer pipeline was then re-run and re-inspected.

Cortical Surface Reconstruction

Cortical reconstruction and volumetric segmentation was performed with the Freesurfer 5.3.0 image analysis suite, which is documented and freely available for download online (http://surfer.nmr.mgh.harvard.edu/). The 3D T1-weighted MPRAGE image was used for the surface reconstruction as it was acquired with the highest spatial resolution. The processing includes removal of non-brain tissue using a hybrid watershed/surface deformation procedure (Ségonne et al. 2004), segmentation of the subcortical white matter and deep gray matter volumetric structures (Fischl et al. 2002, 2004), tessellation of the gray matter white matter boundary, automated topology correction (Fischl et al. 2001; Segonne et al. 2007), and surface deformation following intensity gradients to optimally place the gray/white and gray/cerebrospinal fluid borders at the location where the greatest shift in intensity defines the transition to the other tissue class (Dale and Sereno 1993; Dale et al. 1999; Fischl and Dale 2000). Once the cortical models were complete, inflated cortical surfaces (Fischl, Sereno, Dale 1999; Fischl, Sereno, Tootell et al. 1999) were registered to a spherical atlas which is based on individual cortical folding patterns to match cortical geometry across subjects (Fischl, Sereno, Dale 1999; Fischl, Sereno, Tootell et al. 1999). This method uses both intensity and continuity information from the entire three dimensional MR volume in segmentation and deformation procedures to produce representations of CT, calculated as the closest distance from the gray/white matter boundary to the gray matter/CSF boundary at each vertex on the tessellated surface (Fischl and Dale 2000). The maps are created using spatial intensity gradients across tissue classes and are therefore not simply reliant on absolute signal intensity. The maps produced are not restricted to the voxel resolution of the original data thus are capable of detecting submillimeter differences between groups. Test-retest reliability was validated for Freesurfer morphometric procedures across scanner manufacturers and across field strengths (Han et al. 2006; Reuter et al. 2012). Segmentations were inspected manually for accuracy, and corrected when needed. Errors usually pertained to misclassification of few voxels near membranes or vessels adjacency. Individual cortical and pial surfaces were inspected for errors, and manual correction through control points was performed. The automated Freesurfer pipeline was then re-run for the surfaces that displayed errors and re-inspected. Individual cortical surfaces were then registered to the group template using Freesurfer’s spherical registration.

Analysis of White Matter T1 Maps

For the characterization of myelin in the human brain in vivo, the quantitative T1 mapping by MP2RAGE (Marques et al. 2010) provides a robust measure strongly driven by myelin content of white matter (Lutti et al. 2014; Stüber et al. 2014; Turner 2015). Although qT1 does not reflect uniquely myelin content (Rooney et al. 2007; Fukunaga et al. 2010; Stüber et al. 2014), it can act as a proxy of myelin. T1 relaxation times, although using different sequences from MP2RAGE, have been previously used to obtain white matter microstructural insights in both normal development (Deoni et al. 2012; Yeatman et al. 2014) and pathological alterations of white matter structure, notably in multiple sclerosis (Parry et al. 2002; Vrenken et al. 2006).

The brain masks estimated by Freesurfer were used to skullstrip the T1 maps. Skull-stripped T1 maps have been normalized to the template using the warp fields estimated during the registration to the template. To avoid regions with only partially overlapping white matter tissue between participants, a conservative white matter mask was estimated by registering to the template the individual white matter segmentation masks estimated in Freesurfer. Registration used the aforementioned warp fields and nearest-neighbor interpolation. The resulting masked T1 maps contained only white matter voxels for 100% of the participants, and were thus used for statistical analysis.

Statistical Analysis

Cortical SA was defined as the area of the gray/white matter boundary. SA group-registered maps were smoothed using Gaussian kernel of full-width at half-maximum of 10 mm kernel. Linear effect of age on SA, was estimated through a general linear model for each vertex using Freesurfer’s tool mri_glmfit, including sex as a covariate of no interest. Multiple comparison correction was performed using a False Discovery Rate (FDR) criterion of alpha = 0.01. The vertex-wise analyses were controlled for multiple comparisons using a FDR criterion with α = 0.01. Pearson’s correlation coefficient were calculated from the Z scores for the surviving vertices. Peak Z, P, and ρ values were extracted for each parcel of the Desikan & Killiany atlas (Desikan et al. 2006). Moreover, mean ρ values for each cortical parcel of the Desikan & Killiany atlas were obtained for each hemisphere.

A voxel-wise, nonparametric permutation-based approach was used to test the correlation between the T1 values within the group white matter mask and age, accounting for sex, using FSL’s randomize tool (http://www.fmrib.ox.ac.uk/fsl) with 20 000 permutations. Results were corrected for FDR using α = 0.001. Note that lower T1 values relate to a higher myelination.

Given the high number of voxel-wise comparisons, a cluster-extent threshold criterion was employed. Monte Carlo simulation with 10 000 iterations with a voxel P value of 0.001 using the 3dClustSim program implemented in AFNI (http://afni.nimh.nih.gov/afni/) resulted in an extent threshold of 19 voxels in order to obtain a cluster-level corrected threshold of P < 0.05. Resulting clusters were projected back to the individual space using the inverse warp fields estimated during the registration to the template.

Whole-Brain Probabilistic Tracking

In order to reduce noise in the diffusion-weighted images, a two-step hybrid image restoration procedure was used, consisting of Wiener wavelet filtering followed by speckle reducing anisotropic diffusion (Lohmann et al. 2010). Subsequently, to correct for geometrical distortions caused by magnetic field inhomogeneities, two images acquired with reversed polarity of phase-encoding gradients were used to estimate a nonlinear diffeomorphic warp field (Ruthotto et al. 2012). The corrected images were registered to the individual T1 image using an affine registration in order to correct for motion artifacts. Finally, a diffusion tensor was fitted to the data and the fractional anisotropy (FA) was computed. In order to define an area without significant fiber crossings, subject-specific mask defining the corpus callosum, obtained from the Freesurfer segmentation, was matched to the diffusion-weighted datasets, eroded, and voxels with corresponding FA values >0.7 were selected. The surviving voxels with single fiber populations were used to compute an accurate diffusion-weighted attenuation profile with the MRtrix 3 package (http://www.brain.org.au/software) (Tournier et al. 2012). Spherical deconvolution, using a maximum harmonic order of 8, has been used with this response profile in order to reconstruct the orientation distribution function of voxels within the individual white matter mask derived from the Freesurfer segmentation.

Anatomically-constrained probabilistic white matter tractography was performed using MRtrix within an anatomically-constrained tracking framework (Smith et al. 2012). Seeding at the interface between gray and white matter for unidirectional tracking was performed using the iFOD2 algorithm and backtracking when the streamline failed to enter one of the acceptable ending tissues. A total of 5 million streamline were selected per participant. The white matter clusters with significant correlation between age and T1 relaxation time were used as an inclusion criteria to select the streamlines running through these clusters to the cortex. In order to plot the terminations of these streamlines on the cortical surface, visitation maps of the endpoints were created. The distribution of the streamlines endpoints on the white/gray matter boundary were transformed into surfaces, registered on the common template, averaged across participants and normalized. Averaged normalized values of fiber terminations count for each cortical parcel of the Desikan and Killiany atlas were then extracted and plotted with SA–Age correlation coefficients for display (Fig. 5).

Figure 5.

Spider plot showing the average visitation probability for clusters with significant correlation between age and T1 values of white matter (in green), and complementary for the regionally corresponding correlation between age and cortical surface area (in red).

Statistical Comparison of Surface Area Expansion and Distribution of Fibers Showing Age Effects in T1 Relaxation Times

In order to directly evaluate the relationship between SA expansion with Age and T1 relaxation times effects with Age, a multivariate analysis on SA and Age, accounting for sex was run using the individual visitation maps as per-vertex regressor. Correction for multiple comparison was performed using Monte Carlo simulations as implemented in Freesurfer, using a with a cluster-forming threshold of α = 0.01, testing against an null distribution of maximum cluster size across 10 000 iterations. A mask from the corrected results has been computed at the group level, registered to the single participants as a surface using Freesurfer’s mri_surf2surf tool and transformed into a volume. Fibers whose terminations ended within the significant clusters of this analysis were analyzed through MRTRIX’s tckstats tool in order to extract the average fiber distribution length for each participant, and averaged across the group.

Results

Cortical Surface Area

We first considered developmental effects of age on cortical SA. Figure 1 displays the Z scores of vertices with a significant positive correlation between SA and age (FDR corrected, P < 0.01). The correlation shows a markedly bilateral pattern, encompassing at the lateral aspect the bilateral insula, lateral orbitofrontal cortex, the inferior, middle and superior frontal gyri, inferior parietal lobule, middle and inferior temporal gyri, parahippocampal cortices and fusiform gyri; at the medial aspect the bilateral medial frontal cortices, cingulate gyri, precunei, cunei, and lingual gyri. Table 1 reports the correlation’s local maxima found in each parcel in accordance with the Desikan & Killiany atlas (Desikan et al. 2006). A complementary analysis on CT is added as Supplementary Material.

Figure 1.

Maps (Z scores) showing areas of significant correlation between surface area and age. The values are displayed on the inflated surface. Large regions of cortical surface area, particularly in frontal and temporal lobes, increase with age. Results are FDR corrected (α = 0.01).

Table 1

Peaks of significant correlation (after FDR correction) between age and cortical surface area within each of the parcels described by the Desikan & Killiany atlas. Uncorrected P-values are reported.

Cortical surface area	Left hemisphere		Right hemisphere
Area	Z	ρ	Z	ρ
Banks STS	4.7656	0.3581	5.5471	0.4119
Caudal Anterior Cingulate	5.0156	0.3755	5.7995	0.4288
Caudal Middle Frontal	4.9748	0.3727	4.4672	0.337
Cuneus	—	—	4.5926	0.3459
Entorhinal	4.5324	0.3417	4.6019	0.3466
Frontal Pole	3.9529	0.3001	—	—
Fusiform	5.2326	0.3905	5.0126	0.3753
Inferior Parietal	4.8013	0.3606	4.6677	0.3512
Inferior Temporal	5.6443	0.4184	5.7009	0.4222
Insula	7.4302	0.5326	7.4717	0.5351
Isthmus Cingulate	5.5445	0.4117	—	—
Lateral Occipital	4.4217	0.3338	4.3318	0.3274
Lateral Orbitofrontal	7.3104	0.5253	6.5348	0.4769
Lingual	4.8386	0.3632	5.1858	0.3873
Medial Orbitofrontal	5.2307	0.3903	5.5714	0.4135
Middle Temporal	5.1158	0.3824	5.7653	0.4266
Paracentral	3.7538	0.2857	4.3088	0.3258
Parahippocampal	5.2187	0.3895	4.9615	0.3718
Pars Opercularis	5.6478	0.4187	5.1052	0.3817
Pars Orbitalis	6.6498	0.4842	6.0332	0.4443
Pars Triangularis	4.607	0.3469	6.3896	0.4676
Pericalcarine	4.6687	0.3513	5.1959	0.388
Post-central	5.5098	0.4094	4.3565	0.3292
Posterior Cingulate	4.6937	0.353	4.9308	0.3696
Precentral	5.5675	0.4133	5.0896	0.3806
Precuneus	5.5776	0.414	—	—
Rostral Anterior Cingulate	5.2188	0.3895	6.2735	0.46
Rostral Middle Frontal	5.8888	0.4348	5.0257	0.3762
Superior Frontal	5.2131	0.3891	6.2589	0.4591
Superior Parietal	5.3329	0.3973	4.3063	0.3256
Superior Temporal	4.8019	0.3606	5.1939	0.3878
Supramarginal	5.4787	0.4073	4.2986	0.325
Temporal Pole	4.0262	0.3055	4.9425	0.3704

White Matter Composition

An extensive significant negative correlation between T1 map values and age was found within large parts of the white matter mask. The negative correlation sign was in line with the hypothesis, as during childhood T1 relaxation times of white matter structures decreases (Yeatman et al. 2014). Significance peaks were found bilaterally in the inferior fronto-occipital fascicle, the inferior and superior longitudinal fascicle as well as in the anterior, superior and posterior corona radiata, and in the corpus callosum. Figure 2 represents the T values of the distribution, thresholded at P < 0.001, FDR corrected.

Figure 2.

White matter T-map of voxels displaying a significant negative correlation between T1 values and age. Results are masked using a P value image thresholded using α = 0.001 (FDR corrected). Results are shown on the group template created from the T1 weighted images aligned to the MNI coordinate system.

Combined Analysis of SA and qT1

The core of our analyses focused on the relation between developmental effects as observed for cortical SA and for white matter composition. The streamlines travelling through voxels that displayed a significant Age–T1 relaxation times correlation were selected, and their endpoints were mapped on the cortical surface. The individual endpoints distributions were registered to the template, averaged, and normalize. Figure 3 displays this distribution of the start- and endpoints of the white matter streamlines that contribute to the significant age effects on T1 relaxation times, projected on the common template surface. A comparison with the observed SA–age effects as displayed in Figure 1 reveals a close resemblance between both results. In order to more directly compare age-related correlations in white matter microstructure and myelin content with age effects in SA, both measures are combined in a common map overlapping both effects (Fig. 4A). Large overlaps are evident in the frontal and temporal cortices bilaterally and in the left parietal cortex. Table 2 describes the peak, the mean and the standard deviation of the Z scores resulting from a multivariate analysis of the SA–age relationship using the individual distributions of the significantly myelinating fibers as a per-vertex regressor. Z scores of the analysis are represented on the template surface in Figure 4B. The SA expansion relationship with the underlying fibers composition development is significant within most of the clusters that display an overlap of the two results. A regionally clustered visualization for both results is represented in Figure 5, which shows for each parcel of the Desikan and Killiany atlas the mean ρ values of SA–age correlation together with the averaged probability of endpoints terminations for said parcel. In Figure 6, the average fiber length distribution of the fibers terminating or originating from the clusters depicted in Figure 4B have been plotted. Considering the pediatric population, fibers contributing to the effect appeared to be mostly medium range (mean length = 42.04 mm, standard deviation = 27.13 mm, median = 36.14 mm).

Figure 3.

Normalized average distribution map of the start- and endpoints of connections traveling through the white matter clusters displaying significant negative correlation between T1 values and age. Results are FDR corrected (P < 0.001).

Figure 4.

A: Visualization of significant effects of age on cortical surface area (blue) and on myelin content (T1 values) (yellow) on the same map. A substantial overlap between these two maps is apparent (green). B: Z scores of the per-vertex relationship between Surface Area expansion and the distribution of myelinating fibers’ termination points.

Table 2

Peaks of significant per-vertex multivariate analysis between Surface Area, distribution of terminations of myelinating fibers and age within each of the parcels described by the Desikan & Killiany atlas. Z values (mean, standard deviation and peak) are reported.

	Left hemisphere			Right hemisphere
Parcel	Z Mean	Z StD	Z Max	Z Mean	Z StD	Z Max
Banks STS	1.9676	1.6207	4.4565	2.4885	1.6645	4.479
Caudal Middle Frontal	1.2823	1.7217	4.7847	0.3326	1.0057	4.1045
Cuneus	0.269	0.8469	3.4072	—	—	—
Entorhinal	3.4191	2.4858	8.5852	2.7849	2.0238	5.5118
Fusiform	0.1781	0.7899	4.8665	1.21	1.6686	4.7892
Inferior Parietal	1.3905	1.6839	5.365	1.573	1.7117	4.962
Inferior Temporal	—	—	—	2.3549	1.7867	5.0169
Isthmus Cingulate	0.9041	1.4247	4.0091	—	—	—
Lateral Occipital	—	—	—	0.7825	1.5096	5.0646
Lateral Orbitofrontal	1.6626	1.9538	5.7795	1.0956	1.6709	5.6519
Lingual	0.0011	0.0538	2.6323	0.4752	1.1353	3.8892
Medial Orbitofrontal	1.4921	2.2227	6.3009	0.4232	1.2818	5.6283
Middle Temporal	0.7385	1.4202	4.5188	2.1581	1.7663	5.0302
Parahippocampal	2.3683	2.2302	6.3858	2.1316	1.8525	5.3444
Pars Opercularis	2.227	1.6724	4.9191	0.9821	1.4572	4.0287
Pars Orbitalis	2.698	1.4289	4.9257	1.6567	1.7659	4.5583
Pars Triangularis	2.1829	1.5082	4.0634	4.0482	1.4785	6.2462
Pericalcarine	0.0096	0.1577	2.6454	1.3507	1.5223	3.5758
Post-central	0.8099	1.5032	4.8769	0.2184	0.83	4.6844
Posterior Cingulate	0.3505	0.9495	3.2028	—	—	—
Precentral	0.7712	1.3524	4.3162	0.4571	1.1055	4.327
Precuneus	1.8657	1.7435	5.0404	0.0542	0.3912	3.2446
Rostral Anterior Cingulate	0.4923	1.2843	4.6089	0.0268	0.3151	4.2284
Rostral Middle Frontal	0.8412	1.4362	5.1803	0.6187	1.2452	4.1643
Superior Frontal	0.6457	1.2714	3.927	—	—	—
Superior Parietal	0.3388	0.9772	4.4203	0.3335	0.9788	4.24
Superior Temporal	0.8021	1.518	5.2202	0.736	1.3275	4.0839
Supramarginal	1.3439	1.9212	5.9787	1.4922	1.8041	5.3143
Frontal Pole	0.4883	1.0294	2.739	—	—	—
Temporal Pole	0.0333	0.3026	2.8177	—	—	—
Insula	2.7436	2.1282	7.0343	2.2276	2.083	6.065

Figure 6.

Average streamline length distribution of the fibers originating or terminating from portions of the cortex displaying significant T1/SA relationship with Age (Fig. 3B).

Discussion

The aim of this study was to improve the characterization of gray and white matter maturation and their relationship during development. To do so, we investigated whether cortical morphometric properties and the composition of white matter are linked. We hypothesized that the cortical SA expansion observed during development might be related to the underlying myelination of the connecting cortico-cortical axons. To our knowledge, this is the first study that shows the concurrence between age-related effects in the cortical surface and in the white matter composition at the whole-brain level.

In order to evaluate age effects in the morphometric properties of gray matter, we obtained cortical SA maps from 3- to 7-year-old children. The effect of age on these maps was linearly modeled, while accounting for sex. At the white matter level, in order to observe age effects in the myelination, we analyzed quantitative T1 maps obtained from a MP2RAGE sequence using the same linear model. We then compared the cortical area expansion map with the visitation map of the streamlines running through the clusters that displayed significant decrease T1 values over the analyzed age range.

Our results indicate an extensive increase in cortical SA with the strongest age-related effects in frontal, temporal, and parietal associative cortices and the cingulum. The SA of lingual gyrus, calcarine sulcus, and cuneus also displayed expansion. At the same time, cortical SA of primary motor and sensory cortices remained relatively stable. According to Li et al. (2013), the expansion of somatosensory cortex, superior temporal gyrus and superior parietal lobule peaks within the first year of life. Our results are largely compatible with this view as they show no significant expansion of these areas after age 3 years. Our results also support previous results that indicated an extensive expansion of cortical SA from 4 to 10 years of age encompassing frontal, temporal, and parietal associative cortices (Amlien et al. 2014).

Compared with the expansion of SA, the magnitude of effect of age on CT was less pronounced as shown in our Supplementary Data. CT decrease was localized mainly within the rolandic operculum, the anterior frontal cortex, the inferior parietal lobule, the middle and posterior cingulate gyrus, the middle temporal gyrus and the medial occipital cortex. The regional distribution of these effects is consistent with the results presented for the 4–10 year of age range analyzed by Amlien and colleagues (2014). Compared with previous studies (Gogtay et al. 2004; Shaw et al. 2008; Raznahan et al. 2011), no significant positive correlation between age and CT was observed within our sample. This is in agreement with other recent studies (Mutlu et al. 2013; Nguyen et al. 2013; Amlien et al. 2014; Mills and Tamnes 2014; Wierenga et al. 2014; Croteau-Chonka et al. 2016). Rather, the peak of CT might be reached already at 2 years of age (Li et al. 2013). Notably, we found that the topographic peaks of age effects in SA increase and cortical thinning are not overlapping. This adds to the body of literature suggesting that these two growth patterns might reflect rather independent processes (Amlien et al. 2014; Lyall et al. 2015). The only exceptions to this pattern were found in the lingual gyrus and the calcarine sulcus, which showed significant cortical thinning at the same time as a marked cortical SA increase.

The general pattern of correlations between age and white matter properties was widespread, suggesting ongoing myelination with strongest effects in superior longitudinal fascicle, inferior longitudinal fascicle, inferior fronto-occipital fascicle, and corona radiata. The cortico-spinal tract did not present significant age-related effects. Tracing the connections that contribute to the significant age effects back to the cortex allowed us to obtain a more reliable representation of the cortical regions whose connections display the decrease in T1 values over age. This analysis revealed the strong involvement of frontal, temporal, and parietal associative cortices, together with the medial wall and the precuneus. Interestingly, age does not appear to have a strong effect on white matter connections of visual cortices. However, processes of pruning are still undergoing in this area (Huttenlocher and Dabholkar 1997), which on the other hand might be reflected in cortical volume decrease.

Overlap of Age Effects in Cortical Morphology and White Matter Composition

In order to identify the connection between age effects in gray matter and white matter maturation, we inspected the overlay between both correlation profiles. Projection of the tract endpoints onto the cortical surface yielded visitation map results that strikingly resembled the cortical surface expansion map. In particular, limbic and associative cortices showed an extensive effect of age in both analyses, while primary auditory and somatosensory cortices did not display this effect in either of the two. To further this point, the SA expansion over age was tested using the individual visitation maps as regressor. Given the strong overlap between the patterns expressed by the morphological age effects of the surface and the microstructural age effects of white matter, as shown in Figure 4A, and the highly significant relationship between the two results, as shown in Table 2 and Figure 4B, it is reasonable to suppose that the cortical surface increase is sensitive to the myelination of cortico-cortical fibers projecting to or from the very cortical regions. It is, however, speculative whether this increase in SA is needed to accommodate for the increase in volume of the axonal myelin sheathing. Although this would be in accordance with a stretching model of brain growth (Seldon 2005) other scenarios are possible. Usage-dependent synaptic pruning (Bourgeois et al. 1994), proliferation of myelin into the neuropil (Sowell et al. 2003) and a combination of several of these factors (Amlien et al. 2014) have been suggested to contribute to developmental changes in cortical parameters.

While we found that the peak clusters of age-related effects in CT display a partial overlap with peak age-related effects in white matter microstructure and myelin (see in Supplementary Material, Table 2 and Fig. 2), this interaction does not overlap with the results yielded by the SA–qT1 relationship. Previous studies indicated that myelination and the processes responsible for cortical thinning are mostly independent, as changes in myelin water fraction at the adjacent white matter does not seem to be the driving factor in cortical thinning processes (Croteau-Chonka et al. 2016), which favors the hypothesis that cortical thinning might not be driven by changes in the cortical myelin of the neuropil, as previously suggested (Sowell et al. 2003).

It is worth mentioning that the SA increase cannot uniquely reflect age-related effects in myelin content of the associated white matter fibers, as there are also regions that do not join the large pattern of overlap of the two results, particularly in the medial occipital cortex. This phenomenon might require an alternative explanation. Even though the qT1 does not show an effect, an increase in the volume of gyral white matter cannot be excluded. How the volume and composition of the gyral white matter change during morphological maturation of the cortex is still unknown.

Limitations

It should be noted that the measure of T1 relaxation time is not completely specific for myelin but depends also on other factors such as the iron content and water content of the tissue (Rooney et al. 2007; Fukunaga et al. 2010; Stüber et al. 2014). However, the substantial impact of myelin and more specifically of myelin membrane lipids on T1 relaxation times has been long established (Koenig 1991). Alternative MRI sequences may be used to obtain a similar but complementary measure of myelin content, e.g., myelin water volume fraction (Deoni et al. 2012). But longer acquisition times disqualify such sequences from being used with non-sedated children.

The modeling of age effects on gray and white matter was employed using a linear model. This was chosen because the age range that contributed to our sample was restricted enough to observe only a portion of the nonlinear developmental trajectories observed in studies which observed greater age ranges (Giedd et al. 1999), and within restricted data linear approximation is deemed parsimonious (Bates & Watts 1988). Finally, it should be pointed out again that the current data describe developmental age-related effects using a cross-sectional investigation. Stronger conclusions could be drawn from longitudinal examinations. Future longitudinal studies are invited to complement our results on age-related effects of brain structure and morphology during development.

Click here for additional data file.