Maturational trajectories of pericortical contrast in typical brain development.

In the last few years, a significant amount of work has aimed to characterize maturational trajectories of cortical development. The role of pericortical microstructure putatively characterized as the gray-white matter contrast (GWC) at the pericortical gray-white matter boundary and its relationship to more traditional morphological measures of cortical morphometry has emerged as a means to examine finer grained neuroanatomical underpinnings of cortical changes. In this work, we characterize the GWC developmental trajectories in a representative sample (n = 394) of children and adolescents (~4 to ~22 years of age), with repeated scans (1-3 scans per subject, total scans n = 819). We tested whether linear, quadratic, or cubic trajectories of contrast development best described changes in GWC. A best-fit model was identified vertex-wise across the whole cortex via the Akaike Information Criterion (AIC). GWC across nearly the whole brain was found to significantly change with age. Cubic trajectories were likeliest for 63% of vertices, quadratic trajectories were likeliest for 20% of vertices, and linear trajectories were likeliest for 16% of vertices. A main effect of sex was observed in some regions, where males had a higher GWC than females. However, no sex by age interactions were found on GWC. In summary, our results suggest a progressive decrease in GWC at the pericortical boundary throughout childhood and adolescence. This work contributes to efforts seeking to characterize typical, healthy brain development and, by extension, can help elucidate aberrant developmental trajectories.

Introduction

The human brain is a remarkably complex organ, with a dynamic development profile. By the fourth week of gestation, the brain’s major cortical structures are formed (Muñoz-Moreno et al., 2016), though it continues to develop and organize new neuroanatomical pathways well into adulthood (Giedd and Rapoport, 2010; Lebel and Beaulieu, 2011; Shaw et al., 2008). Cortical maturation and the refinement of neuronal pathways have been shown to have direct links to children’s behavior and cognitive performance and potentially contribute to identify maladaptive trajectories of development and associated neuropsychiatric disorders (Ducharme et al., 2014; Ducharme et al., 2012; Giedd et al., 2009; Giedd and Rapoport, 2010; Kharitonova et al., 2013; Lenroot et al., 2007; Shaw et al., 2011; Shaw et al., 2006, 2006; Thormodsen et al., 2013). Improving our knowledge of the neural development of healthy children and adolescents may hence allow us to identify atypical trajectories as well as critical developmental periods when a child’s brain is most vulnerable to environmental influences.

Over the last three decades, neuroimaging has become an indispensable tool in brain research, allowing for in-vivo, non-invasive characterization of neurological structures at a submillimeter level of resolution (Ameis et al., 2014; Dosenbach et al., 2010; Giedd et al., 1999; Gogtay et al., 2004; Khundrakpam et al., 2016; Raznahan et al., 2011; Shaw et al., 2006, 2006; Thompson et al., 2000). Significant efforts have been made to accurately describe the developmental trajectory of several morphological features of the cortex as measured with magnetic resonance imaging (MRI). Cortical volume was one of the first cortical metrics to be longitudinally characterized and was reported, in earlier studies, to increase until late childhood and then to decrease afterwards (Giedd et al., 1999; Gogtay et al., 2004; Lenroot et al., 2007). Later studies have reported that cortical volume reached its peak in early childhood and decreased in later childhood and throughout adolescence (Lebel and Beaulieu, 2011; Mills et al., 2016; Mills et al., 2014; Tamnes et al., 2017, 2013).

In order to more precisely characterize cortical development, cortical volume measures can be effectively fractionated into cortical thickness and surface area components, which have been shown to be more susceptible to environment and genetics, respectively (Ducharme et al., 2015; Fischl et al., 1999; Ghosh et al., 2010; Patel et al., 2020; Postelnicu et al., 2009; Raznahan et al., 2011; Thompson et al., 2017). There have been inconsistencies between studies regarding the trajectory of cortical thickness development from mid-childhood to late adolescence (Tamnes et al., 2017). Some early work (Raznahan et al., 2011; Shaw et al., 2007) reported increasing thickness until late childhood followed by a decline afterwards while the majority of studies reported a general decrease with age during mid-childhood and adolescence (Alexander-Bloch et al., 2014; Ducharme et al., 2016; Fjell et al., 2015; Mills et al., 2014; Sowell et al., 2004; Tamnes et al., 2017; Zhou et al., 2015; Zielinski et al., 2014). As for surface area, the majority of studies report a tendency for an increase during childhood followed by a decrease during adolescence (Ducharme et al., 2015; Mills et al., 2014; Raznahan et al., 2011; Tamnes et al., 2017).

Recently, methods to measure T1-weighted (T1w) pericortical gray-to-white matter intensity ratio (referred to as “gray-white contrast”/GWC in this work), have been developed; GWC can be calculated as the ratio of the T1w signal intensity at a certain distance on either side of the gray-white cortical boundary. Some have used the gray to white ratio as their measure of contrast (Salat et al., 2009), with larger values representing less sharp pericortical boundaries, whereas others have instead used white to gray ratio as their measure of contrast, with larger values representing sharper pericortical boundaries (Lewis et al., 2018). In the current work, we adopt the white to gray ratio approach, with larger values representing sharper boundaries.

GWC is considered to be influenced by the development of pericortical microstructure including myelin invading the lower layers of the cortex during development (Dale et al., 1999; Grydeland et al., 2013; Norbom et al., 2019; Patel et al., 2020; Rowley et al., 2015; Vidal-Piñeiro et al., 2016). As such, GWC is thought to provide insight into subtle variations of neuronal tissue that are not otherwise captured by cortical thickness, area, or volume metrics alone. There are nonetheless likely certain relationships between GWC and these other MRI-based metrics. Indeed, the placement of the gray-white surface boundary is based on GWC and may hence influence estimates of cortical thickness, area, and volume. In fact, the invasion of myelin in the lower layers of the cortex may be responsible for part of the apparent reductions in MRI-based estimates of cortical thickness previously attributed exclusively to gray-matter processes. This sentiment is also put forth by Sowell et al. (2004), wherein they note that the invasion of the cortex by myelin may yield lower estimates of MRI-based cortical thickness.

Other studies that have worked on discerning the contributions of biological processes which drive maturational changes in the cortex have considered microstructural pericortical measures (such as GWC) as distinct and informative measures (Croteau-Chonka et al., 2016; Deoni et al., 2015; Norbom et al., 2019, 2020). For example, Croteau-Chonka et al. (2016) found that maturation of cortical thickness in early childhood in some regions of the cortex was unaffected by white-matter maturation, concluding that brain imaging studies of cortical and white matter maturation can reflect distinct, but complimentary, neurodevelopmental processes. Importantly, GWC has previously been established as a metric linked to cognitive performance (Lewis et al., 2018), aging of the brain (Magnaldi et al., 1993; Raz et al., 2005), and developmental neuropsychopathology (Bezgin et al., 2018). Still, there is a relative paucity of studies in this domain, especially those seeking to disentangle the relative contributions of biological processes on either side of the pericortical gray-white boundary. The current study follows directly from previous efforts examining longitudinal pericortical maturation by providing a detailed characterization of the developmental trajectories of GWC in healthy children and adolescents (Brain Development Cooperative, 2012). This is the first study that uses repeated scans within a longitudinal design to examine GWC changes across the important maturational changes that occur in typically developing individuals between childhood and early adulthood.

Materials & Methods

Participants

Data were obtained from the longitudinal NIH Study of Normal Brain Development (Evans and Brain Development Cooperative, 2006). Study participants were recruited across six Pediatric Study Centers across the United States and ranged in age from 4.5 to 18.3 years (N = 431 participants; 207 M/224 F) at the first visit. They were selected using a population-based sampling plan aiming to recruit a sample representative of the US population, based on data from the 2000 US census (Brain Development Cooperative, 2012). More specifically, the sampling plan allowed to build a sample representative of the general US population (2000 census) by referencing geocoded census data. More specifically, recruitment was monitored continuously to maintain a demographically representative sample based on age, gender, ethnicity, and socioeconomic status variables (Evans and Brain Development Cooperative, 2006). Subjects were interviewed following the Diagnostic Interview for Children and Adolescents. Those found to have a major Axis I diagnosis according to DSM-IV classification, were excluded. Further, children with a Child Behavior Checklist syndrome T-score equal to or greater than 70, abnormal findings on neurological examination, or exposure during gestation to substances that may alter brain development, were also excluded. The full list of selection criteria can be accessed in Evans and Brain Development Cooperative (2006). Subjects underwent repeated magnetic resonance brain imaging (MRI) every 2 years, with a maximum of three scans over a 4-year period.

MRI protocol

Whole-brain 3D T1-weighted spoiled gradient recalled echo sequences were acquired on 1.5T MRI scanners across sites with a 1 mm isotropic resolution except for GE scanners which had an in-plane 1 mm isotropic but for which a 1.5 mm slice thickness was allowed due to their limit of 124 slices (Evans and Brain Development Cooperative, 2006). Inter-site reliability of anatomical measurements was assessed throughout the study by scanning phantoms of the American College of Radiology, as well as live volunteers regularly and at each site (Evans and Brain Development Cooperative, 2006; Kabani et al., 2001; Lerch and Evans, 2005). Scans were reviewed immediately after acquisition and repeated if there were either significant motion artifacts or field distortions. As described in Evans and Brain Development Cooperative (2006), only scans that passed raw MRI quality control were included in the final release of the NIHPD dataset. Another quality control (QC) procedure was implemented after image processing (see Section 2.4 below).

Image processing

The CIVET 2.1.1 pipeline developed at the Montreal Neurological Institute was used to extract cortical surfaces on T1w images for corticometric analyses (http://www.bic.mni.mcgill.ca/ServicesSoftware/CIVET-2-1-0-References). In order to take into consideration brain volume variations between subjects, the native MR images were linearly registered to standardized MNI-Talairach space, based on the MNI-ICBM152 non-linear symmetric model (Collins et al., 1994; Fonov et al., 2009; Grabner et al., 2006; Mazziotta et al., 1995; Talairach and Tournoux, 1988). Intensity inhomogeneity introduced by the scanner was normalized via N3 (Sled et al., 1998). The images were then classified into cerebrospinal fluid, gray matter, white matter, and background image by a neural net classifier (INSECT) (Tohka et al., 2004; Zijdenbos et al., 2002). Cortical surfaces were derived with a marching-cubes algorithm (Lorensen et al.), ultimately demarcating 81,924 vertices on every subject’s brain. Vertex correspondence is done via a Laplacian map method described in Kim et al. (2005). In the subjects’ native space, new gray and white matter sampling surfaces were generated using the CIVET/2.1.1 function average_surfaces, which creates a surface whose vertex coordinates are the average of the vertex coordinates of the two input surfaces. The surface at 25% of the distance from the gray/white boundary to the pial boundary (‘25% surface’) was generated by first averaging the gray/white boundary and pial surfaces to create a surface that is 50% of the cortical distance, and then averaging the gray/white surface with the 50% surface. New white matter surfaces were generated at the same distance from the gray/white boundary surface as the gray matter surfaces, but in the direction of the white matter. Specifically, each new white matter surface was generated by first creating vectors between each vertex on the gray/white boundary surface and the corresponding vertex on the gray matter surfaces (separately for each gray matter surface generated). The vectors were then inverted in the x, y, and z dimensions and added to the coordinates of the gray/white boundary surface vertices, such that white matter surfaces were created at the same distance from the gray/white boundary as the gray matter surfaces, but in the opposite direction (toward the superficial white matter instead of toward the pial surface). To establish correspondence of these vertices between subjects, the cortical surfaces were registered non-linearly to a high-resolution average surface template of the ICBM152 data set (Boucher et al., 2009; Lyttelton et al., 2007; Robbins et al., 2004).

Quality control of processed images

The sensitivity of automated cortical measurements to scan-specific field distortions and to small movement artifacts has already been established (Reuter et al., 2015), and the need for a quality control of pipeline outputs is now considered necessary as failing to do so can significantly alter findings when it comes to developmental trajectories (Ducharme et al., 2016; Olafson et al., 2020).

As a first step, we applied the automated CIVET QC protocol which produces, for each subject, multiple quantitative metrics including number of self-intersecting and gray-white surface intersections, brain mask extraction errors, and percentage of error in expansion of gray matter, among others: http://www.bic.mni.mcgill.ca/ServicesSoftware/CIVET-1-1-12-Quality-Control. After excluding problematic outputs, we then applied the same quality control procedure used by Ducharme et al. (2016). Two investigators (ACT and SK), blinded to all clinical variables, implemented quality control of the pipeline outputs by means of unanimous consent, according to previously established quality control criteria; note that SK was a co-investigator in the Ducharme et al. manuscript (2016). Specifically, processed scans were removed if they featured gross deformities, cropped areas in surfaces, and/or poor delineation between gray and white matter surfaces in large areas; for an example comparison between failed scans and an accepted one, see Supplementary Fig. 1.

In the process of QC, 113 scans were removed. From the remaining scans, some had missing information at visits, such as scanner identification or age, which amounted to another 90 scans to be removed. Additionally, prior to analysis, individual scans were removed if whole-brain mean contrast was more than three standard deviations from the mean of the sample - this led to removing only two scans, each from a different subject. A total of 819 scans from 394 participants were retained. There were 394 subjects with one timepoint, 290 subjects with two timepoints, and 135 subjects with three timepoints. Supplementary Fig. 2 shows the age distribution of our dataset. Under the NIHPD terminology, 296 participants were scanned at Visit 1, 313 participants were scanned at Visit 2, and 210 participants were scanned at Visit 3.

Gray-White contrast

GWC is defined here as the ratio of the intensities of white and gray matter on T1w MRI images after N3 correction. In the subjects’ native space, we created a gray surface at 25% of the distance from the gray/white boundary to the pial boundary and a white surface at the same distance but in the direction of white matter. A relative distance to compute GWC was chosen in favor of a fixed distance due to the presence of extremely thin sections of cortex, wherein a fixed distance could, in some instances, lead to sampling at a distance greater than the thickness of the cortex. A relative sampling distance was also preferred over a fixed one to maximizes chances of sampling from the vicinity of the same cortical layer(s) across the cortical mantle. The value of 25% was based on results from Whitaker et al. (2016), showing that the highest myelin concentration observable in the cortex is within the 20% to 30% range of local cortical thickness. To produce the final contrast measures, the sub-white surface values were divided by the gray surface values at each vertex; thus, a larger GWC value represents a sharper gray-white boundary. An example depicting the placement of the surfaces used to compute contrast at each vertex is found in Fig. 1. Finally, spatial surface blurring was applied using a Gaussian blur with FWHM of 20 mm.

Fig. 1.

Calculation of GWC, Example showing the placement of the 25% distance surfaces in the gray matter and superficial white matter.

Assessment of partial volume effects

Given that a voxel may be composed of a mixture of tissue types, an effect known as partial volume effect (PVE), and that such issues may have corollary consequences on our estimations of contrast, we examined PVE on our data. To do so, we used the trimmed minimum covariance determinant (TMCD) method as described by Tohka et al. (2004) for the estimation of the commonly used mixel PVE model (Choi et al., 1991). TMCD provides accurate and robust estimation of partial volume parameters that can be used for the correct delineation of cortical surfaces (Tohka et al., 2004). Both the gray and white surfaces were examined for PVEs from surrounding tissue.

Statistical analyses

Whole-brain, vertex-wise analyses were conducted using Surfstat (https://www.math.mcgill.ca/keith/surfstat/), a statistical toolbox implemented in MATLAB. Mixed-effects model (with intercept as the random effect for each subject) were fit using residualized maximum likelihood estimation (REML). These mixed-effects regression models, which control for within-subject variance, can deal with multiple measures per subject, missing data, as well as irregular intervals between visits. By adding a random intercept in the model, we implicitly assume that observations within subjects are correlated, with a fixed correlation for each pair of observations. This model also allows to model subject-specific age trajectories, with a constant slope (for linear, quadratic, and cubic age effects) across subjects and subject-specific intercepts. Note that we don’t have enough data points over time to include a random slope (and therefore a subject-specific slope) for the age coefficient. Non vertex-wise analyses were performed in R version 3.6.2 using the nlme package (Pinheiro et al., 2013). Results were thresholded using a 5% false discovery rate (FDR).

The following three models were applied across 77,122 vertices (brainstem and non-cortical midline are masked out) and the Akaike Information Criterion (AIC) was used to select the best fitting model for each vertex (Cavanaugh and Neath, 2019).

The inclusion of the age variable allows for the implicit assessment of rate of change with time. The potential for interactions between age and sex was examined using SurfStat. Mixed-effects modeling was performed as described above, except with the inclusion of the appropriate Sex:Agen interaction terms, again, based on the lowest AIC. Sex* Age, Sex* Age2, and Sex* Age3 interactions were examined within full models including lower order interaction(s) for Sex* Age2 and Sex* Age3. To showcase GWC trajectories in a variety of regions across the brain, mean GWC was calculated within a 3 mm radius circle along the surface (mean included vertices = 15) within each selected region. A radius of 3 mm was chosen in order to sample from a large enough area while avoiding sampling from different trajectory types (i.e. avoid examining a region where the center of the 3mm-radius circle is within a region with a quadratic trajectory type while the perimeter is within a region with a linear trajectory type). To communicate regionality, the Desikan surface atlas labels were used (Desikan et al., 2006; Klein and Tourville, 2012). Mixed-effects models were fit as described above, with mean GWC as the dependent variable; the model type was determined by the prevailing model in the selected region, from the aforementioned whole-brain, vertex-wise analysis. To produce the scatterplots of each locale’s progression, mean GWC was predicted from the model fits; modeling was performed independently for male- and female-only subsets of the data, in addition to the global curve derived from the full sample.

Finally, to visualize the dynamics of how maturation of GWC progresses across the brain, vertex-wise estimations of GWC trajectory were derived and mapped onto the cortical surface. Specifically, mixed-effects models were fit at each vertex, again based on lowest AIC from the original vertex-wise analysis. Two methods were used to help visualize, through time, the progression of GWC; in the first, the predicted GWC values were mean centered, and in the other, they were standardized via vertex-wise scaling between 0 and 1 using the R function scales::rescale. These modified vertex-wise GWC predictions were then color-coded and mapped onto the brain surface. A slight 3 mm smoothing was used to improve visualization. Five snapshots were chosen to depict progression throughout childhood and into late adolescence.

Results

Demographics

No statistically significant differences were found between the original sample and our quality-controlled subsamples in terms of gender distribution, ethnic representation and socioeconomic status (see Table 1 and Supplementary Tables 1 and 2). However, our Visit 1 subsample had a slightly higher mean age than the original NIH Visit 1 sample (t = 3.17, p ≤ 0.002). This is because younger children had slightly higher QC failures likely due to having a more difficult time remaining still in the scanner (Reuter et al., 2015).

Table 1

Sample comparison between the original NIHPD Objective 1, Visit 1 sample, and the Visit 1 sample used in this study. When appropriate, means ± standard deviations are provided.

Sample	Objective 1 Visit 1 (n = 433)	Accepted Visit 1 Sample (n = 296)	Statistics (if applicable)
Age (Years)	10.4 ± 3.8 (4.6 to 18.3 years)	11.3 ± 3.7 (4.9 to 18.3 years)	t = 3.17, p < 0.002
Proportion of Males	48%	45.2%	χ2 = 0.301, p = 0.583
Proportion of low/medium/high adjusted SES[a]	24.4% / 40.8% / 34.8%	22.6%/41.2%/36.1%	χ2 = 0.185, p = 0.912
Proportion of Caucasians/African American/Other[b]	73.5% / 8.4% / 18.1%	74.7% / 8.4% / 16.9%	χ2 = 0.099, p = 0.951

Based on the US Department of Housing and Urban Development method for comparing family income levels as a function of regional costs of living.

The ‘Other’ category includes American Indian, Alaskan Native, Asian, Native Hawaiian or Other Pacific Islander, and those for whom ethnicity or race was not provided, or whose parents came from different racial or ethnic backgrounds.

Results of the examination of partial volume effects

With a few exceptions, there was little evidence for significant PVE impacts on our results. This is seen in Supplementary Fig. 3, where average PVE estimates from surrounding tissue are mapped onto the gray and white surfaces that were used to estimate GWC. Not finding strong evidence for PVEs is compatible with the fact that the lowest mean thickness across the cortex, for the sample of scans in this work, was above 1.5 mm (i.e., above the 1 mm isotropic voxel size used in the vast majority of scanners in NIH Study of Normal Brain Development study). For the gray surface, essentially no cerebrospinal fluid PVE was detected. However, a few regions exhibited some white matter PVE effects. These included the dorsal pericentral areas as well as the medial and posterior lateral occipital regions, bilaterally. This was to be expected because these regions are known to exhibit deep myelin penetration (Carey et al., 2018). As for the white surface, some PVE was detected in the insular region, bilaterally.

Overview of gwc findings

As seen in Fig. 2, GWC undergoes widespread developmental change. Vertex-wise, mixed-effects modeling produced a map wherein approximately 99% of vertices attained significance under linear, quadratic, or cubic trajectories of development. Approximately 16% of vertices reflected linear developmental trajectories of GWC as likeliest, 20% depicted quadratic trajectories as likeliest, while 63% of vertices reflected cubic trajectories as likeliest. Overall, most trajectories evidenced trends of declining GWC with age. Models that showed trajectories of increasing GWC were sparse, with slopes close to zero, and were exclusively located in the temporal poles and the medial orbitofrontal cortex. The quadratic developmental trajectories all appeared to be concave downward. Trajectories were also remarkably symmetric between hemispheres, with a few exceptions, outlined below.

Fig. 2.

Map of vertex-wise Model Winners of GWC trajectories, 5% FDR.

Prevailing model based on AIC denoted: Blue = Linear, Green = Quadratic, Purple = Cubic.

GWC findings on whole sample

Whole-brain mean GWC was found to follow a cubic developmental trajectory based upon the lowest AIC (see Fig. 3); this trajectory was shown to follow a downward/decreasing trend with time. Our vertex-wise estimates parallel this finding with the main GWC trajectory being cubic. The distribution of the likeliest trajectories (i.e. linear, quadratic or cubic) were highly symmetrical and showed negative associations between GWC and age across the majority of the brain (See Fig. 2 and Supplementary Fig. 4).

Fig. 3.

Cubic age trajectory of whole-brain, mean GWC.

Black curve is derived from the model prediction on the full sample. The blue and red curves represent predictions from modeling performed on the male- and female-only subsets of the data, respectively.

Modeling the estimated, region-specific developmental progression of GWC in a handful of representative regions (Figs. 4a and 4b) yielded several findings. Qualitative assessments of trajectories within regions depicting quadratic and cubic trajectories of development showed many areas where GWC declined at an increasing rate between the ages of 10 and 15 years and continued its decline well into late adolescence. The temporal, frontal and parietal lobes all mostly exhibited cubic trajectories of decreasing GWC, while the occipital lobe evidenced a mixture of linear and cubic downward trajectories. In the limbic regions, developmental trajectories were mixed between cubic and some linear trends. An example of this is the cingulate, with the anterior cingulate reflecting a linear trajectory, while the posterior cingulate was found to follow a cubic trend. The primary motor cortex also included some linear declining trajectories. Notably, the cuneus exhibited a mix of linear and cubic decline. The superior parietal lobule, posterior temporal lobe, and all the precuneus, bilaterally followed linear trajectories of decline. There were also two symmetrical, anterior-posterior strips on the superior frontal gyri, which were found to deviate from the predominant cubic trajectory in the surrounding region of the frontal lobe; the strip in the left superior frontal gyrus evidenced a linear maturational trajectory, while the strip in the right superior frontal gyrus was observed as having a quadratic trajectory.

Fig. 4.

A. Qualitative overview of developmental trajectories of GWC, on the left lateral cortical surface. Black curves are derived from the model prediction on the full sample. The blue and red curves represent predictions from modeling performed on male- and female-only subsets of the data, respectively, and better reflect the true trajectories of males and females in our data. Even if these are not perfectly parallel to the global trajectory, there is no statistically significant ‘sex by age’ interaction. B. Qualitative overview of developmental trajectories of GWC, on the left medial cortical surface. Black curves are derived from the model prediction on the full sample. The blue and red curves represent predictions from modeling performed on male- and female-only subsets of the data, respectively, and better reflect the true trajectories of males and females in our data. Even if these are not perfectly parallel to the global trajectory, there is no statistically significant ‘sex by age’ interaction.

While the majority of the brain was found to follow a trajectory of decline in GWC, a few noteworthy exceptions are reported as follows. The right lateral occipital lobe depicts increasing GWC, as does the left posterior cingulate. Additionally, there appears to be a slight delay in progression of GWC to its minimum in the orbitofrontal cortex, relative to that of the frontal cortex as a whole. Overall, the frontal cortex was estimated to reach its minimum GWC later (within our age range), when compared to the rest of the brain (See Supplementary Figures 5 and 6).

Sex differences

Vertex-wise analyses that examined the effect of sex on GWC yielded a few significant regions where females tended to have lower contrast values than males across the age range examined here. Regions with a significant sex contrast (males-females) included the cuneus and parts of the precuneus, the anterior cingulate, and a large section of the medial orbitofrontal cortex. Fig. 5 shows a selection of plots produced from the significant regions on the left hemisphere (see Supplementary Figure 7 for areas of significant sex effects on both hemispheres); estimated trajectories show females consistently exhibiting a lower GWC. No interaction between age and sex was found to be statistically significant.

Fig. 5.

Estimated trajectories of GWC maturation across various representative regions exhibiting a significant sex effect on GWC.

Black curves are derived from the model prediction on the full sample. The blue and red curves represent predictions from modeling performed on male- and female-only subsets of the data, respectively, and better reflect the true trajectories of males and females in our data. Even if these are not perfectly parallel to the global trajectory, there is no statistically significant ‘sex by age’ interaction.

Discussion

This work constitutes the first longitudinal study of GWC neurotypical development in a representative American population during the critical developmental period between childhood and early adulthood (Waber et al., 2007). The vertex-wise analysis of the changes in GWC during the development of children and adolescents in our sample yielded a mixture of linear, quadratic, or cubic trajectories as significant across nearly the whole cortical surface. The cubic trajectory was the one that was most frequently observed. Overall, these results predominantly showed globally declining GWC with age. We also observed that the whole-brain average GWC followed a declining cubic trajectory. This extends prior findings from Vidal-Piñeiro et al. (2016), who found GWC to decline from early adulthood to old age. Our findings also fit well with prior work from (Norbom et al., 2019), who found GWC to globally decrease with age in a cross-sectional analysis of a sample of subjects that included various types of psychopathology. When viewing the vertex-wise, regional, and average whole-brain GWC trajectories, the complexity of maturational changes in GWC throughout early-life development is apparent. A prominent finding was the remarkable symmetry across much of the cortical surface, regarding the best fitting model as seen in Fig. 2 and Supplementary Figure 4.

GWC in the frontal cortex appeared to show a delayed maturation as it manifested a decrease in contrast later than in the other regions, possibly reflecting the prolonged penetration of intracortical myelin (relative to the rest of the cortex) known to take place well past our examined age-range in this region (Bartzokis, 2004; Grydeland et al., 2019; Westlye et al., 2010). Speculation may suggest that healthy maturation of certain regions requires maintaining sharp boundaries between gray and white matter (i.e., high GWC) well into adulthood and old age. This view is compatible with results from Vidal-Piñeiro et al. (2016), who showed a positive association between symptoms of aging and lower GWC in the frontal cortex. Alternatively, Vidal-Piñeiro et al’s (2016) results could simply reflect a higher starting point of GWC in early life (followed by a decrease in adulthood) for those having greater cognitive reserve.

In the cingulate cortex, our results quite elegantly show an obvious demarcation between the anterior and posterior cingulate, which fits well with the known cytoarchitectural and functional (e.g. error detection and emotional processing, respectively) roles of these two cingulate regions (Lane et al., 1998; Leech and Sharp, 2014; Leech and Smallwood, 2019; Stevens et al., 2011). A similar finding is visible on the insula (See Supplementary Figure 4), where the inflated cortical surface clearly shows a demarcation that fits well with the known histological demarcation between the anterior agranular insula and the more posterior granular insula (Afif et al., 2010; Naqvi et al., 2014).

Our analyses evaluated the effect of sex on GWC, as well as whether there was a significant interaction between age and sex on GWC. Sex was found to have a significant effect (at 5% FDR) on GWC in several brain regions, across all three types of age trajectories. As shown in Fig. 5, males presented with slightly higher GWC than females in regions where a significant main effect of sex was found. No significant interaction effect on GWC was found anywhere across the cortical mantle for Age* Sex, Age2*Sex or Age3*Sex interactions. In keeping with this, male and female trajectories showed strikingly similar trajectories even though the plots were produced independently for the male and female samples.

The main effect of sex that is observed could potentially be linked to previously reported differences between males and females in both gray and white matter as informed by prior works. For instance, a study from Raznahan and colleagues found a higher percentage of gray matter tissue and interhemispheric connectivity in females, and higher percentage of white matter tissue and intrahemispheric connectivity in males (Raznahan et al., 2011). Ingalhalikar et al. (2014) put forth strong evidence of differences in structural connectivity between males and females, with their results indicating that male brains are structured to facilitate intrahemispheric connectivity when compared to females, while female brains conversely had higher interhemispheric connectivity and cross-modular connections that would facilitate interhemispheric connectivity. Considering this, we propose that the widespread effect of sex on GWC reflects sex differences in intra-cortical patterns of myelination that support sex differences in connectivity profiles. In other words, we speculate that short- (intrahemispheric) and long- (interhemispheric) range connectivity is reflected by myelination patterns that lead to high and low GWC, respectively.

While the overall trend is of decreasing GWC, parts of the right lateral occipital lobe, temporal poles, parahippocampal gyrus, and orbitofrontal cortex were exceptions and exhibited trajectories of increasing or essentially level GWC. Given potentially altered MRI signal due to magnetic susceptibility between air-filled sinuses in the vicinity of some of these brain regions (Juchem et al., 2010) and/or to possible partial volume effects, both of which may have been too subtle to detect by QC in some cases, we consider the isolated findings of increasing GWC in some of these regions to potentially be artefactual (Du et al., 2007; Kharabian Masouleh et al., 2020; Reuter et al., 2015; Savalia et al., 2017). Regarding partial volume effects, it is noteworthy that such effects were shown to be present in the insula, uncus, parahippocampal gyrus (as mentioned above), dorsal pericentral, medial occipital, and posterior lateral occipital regions and making estimates of GWC in these regions as tentative.

It is important to keep in mind that the biological underpinnings of GWC and of variations in GWC are not yet clearly elucidated (Lewis et al., 2018; Raznahan et al., 2011; Salat et al., 2009). As our dataset lacks true data representing myelinated axons at the cortical boundary, we are bound to speculate on underpinnings of the negative trajectories of GWC change observed here, while acknowledging the partial contribution from other sources that would affect GWC makeup, such as cellular iron content (Fukunaga et al., 2010). One possibility is that GWC changes reflect the selective pruning of gray matter in such a way that it would affect the ratio value of gray to white matter intensity (Elston et al., 2009). Another, perhaps more plausible alternative that is gaining traction is that contrast changes rather reflect the gradual invasion of the lower layers of cortex by myelin (Elston et al., 2009; Kharitonova et al., 2013; Petanjek et al., 2011; Rakic, 2002). If this is the case, our data may very well reflect genuine myelination of gray-matter fibers in the lower layers of the cortex. This assertion is supported by closer examinations of early life development around our given age-range (Deoni et al., 2015; Gilmore et al., 2012; Sowell et al., 2003). Indeed, various neuroimaging studies that used volumetric methods to measure growth of gray and white matter have shown a gradual decrease in gray matter volume that coincides with an increase in white matter volume (Courchesne et al., 2000; Giedd et al., 1999; Jernigan et al., 1991; Sowell et al., 2003; Sowell et al., 2004; Sowell et al., 2001). This is further supported by postmortem studies that show myelin penetration well into adolescence and a lack of significant loss of cortical neuronal cell bodies until the fourth or fifth decade of life (Dekaban and Sadowsky, 1978; Raznahan et al., 2011; Terry et al., 1987; Wang and Young, 2014), as well as by histological data that shows the cortical boundary as anything but discretely “black or white” (Annese et al., 2004). In fact, what is currently referred to as “cortical thinning” in early life may largely be attributed to myelination of neuronal fibers in the lower layers of the cortex, as opposed to physically smaller or fewer gray-matter neurons and associated processes (Sowell et al., 2004). In sum, myelination of intracortical gray matter fibers is likely responsible for most of the apparent reduction in GWC found in our data.

Having said this, we must also acknowledge that GWC could be calculated differently. For example, while our measure of GWC was based on a relative distance from the gray/white matter boundary (e.g., 25% distance from gray-white surface boundary to the pial surface), a fixed distance (e.g., 1 mm) or another relative distance could have been chosen instead. Alternatively, a relative distance within the gray-matter, and a fixed distance (e.g., 1 mm) into the white matter could have been possible as well. As discussed in the methods, for the gray surface, we favored a sampling distance that was relative to gray-matter thickness (instead of a fixed distance) to decrease the possibility of sampling from regions outside the cortex in areas of very thin cortex.

A caveat of the current work is that our age range may have influenced our observations. We propose, for instance, that quadratic or cubic trajectories of decline would emerge as likeliest in most regions where only a linear trend was observed, because the slope of a linear decrease is likely to attenuate at some point with age to avoid reaching a GWC of zero early in adulthood. This reasoning can also be extended to regions of quadratic trajectories that might follow a cubic trajectory had a larger age range been available in our cohort. Another caveat pertains to partial volume effects.

Ultimately, additional, and more expansive investigations will be required for elucidating the nature of GWC changes across the lifespan. For instance, future work on GWC would benefit from using scans with sub-millimeter voxel sizes readily achievable with higher field strength (e.g., 3T). Such work would decrease potential partial volume effects, improve surface extraction, and could help improve our understanding of how GWC combines with other cortical metrics (e.g., cortical thickness). Also, further refinements in how GWC is calculated and using GWC to adjust estimates of other cortical metrics might help improve our characterization of brain development.

In summary, we find that the distribution of the best fitting GWC trajectory model is highly symmetrical across hemispheres and that the GWC tends to decrease with age during development over most of the brain. While the underlying neurobiology of GWC is not entirely clear, the current work, which constitutes the first longitudinal GWC study of the critical developmental period between childhood and early adulthood, adds to the growing body of knowledge on cortical development (Bartzokis, 2004; Deoni et al., 2015; Glasser et al., 2014; Glasser and Van Essen, 2011).