Abstract
Adolescence is characterized by the maturation of cortical microstructure and connectivity supporting complex cognition and behavior. Axonal myelination influences brain connectivity during development by enhancing neural signaling speed and inhibiting plasticity. However, the maturational timing of cortical myelination during human adolescence remains poorly understood. Here, we take advantage of recent advances in high-resolution cortical T1w/T2w mapping methods, including principled correction of B1+ transmit field effects, using data from the Human Connectome Project in Development (HCP-D; N = 628, ages 8–21). We characterize microstructural changes relevant to myelination by estimating age-related differences in T1w/T2w throughout the cerebral neocortex from childhood to early adulthood. We apply Bayesian spline models and clustering analysis to demonstrate graded variation in age-dependent cortical T1w/T2w differences that are correlated with the sensorimotor-association (S-A) axis of cortical organization reported by others. In sensorimotor areas, T1w/T2w ratio measures start at high levels at early ages, increase at a fast pace, and decelerate at later ages (18–21). In intermediate multimodal areas along the S-A axis, T1w/T2w starts at intermediate levels and increases linearly at an intermediate pace. In transmodal/paralimbic association areas, T1w/T2w starts at low levels and increases linearly at the slowest pace. These data provide evidence for graded variation of the T1w/T2w ratio along the S-A axis that may reflect cortical myelination changes during adolescence underlying the development of complex information processing and psychological functioning. We discuss the implications of these results as well as caveats in interpreting magnetic resonance imaging (MRI)-based estimates of myelination.
SIGNIFICANCE STATEMENT Myelin is a lipid membrane that is essential to healthy brain function. Myelin wraps axons to increase neural signaling speed, enabling complex neuronal functioning underlying learning and cognition. Here, we characterize the developmental timing of myelination across the cerebral cortex during adolescence using a noninvasive proxy measure, T1w/T2w mapping. Our results provide new evidence demonstrating graded variation across the cortex in the timing of T1w/T2w changes during adolescence, with rapid T1w/T2w increases in lower-order sensory areas and gradual T1w/T2w increases in higher-order association areas. This spatial pattern of microstructural brain development closely parallels the sensorimotor-to-association axis of cortical organization and plasticity during ontogeny.
Introduction
Adolescence, the transition from childhood to adulthood, is characterized by the refinement and stabilization of neural circuits supporting the dynamic control of attention and behavior (Larsen and Luna, 2018). Axonal myelination in the cerebral cortex enhances signaling speed (McDougall et al., 2018) and regulates synapse formation during development to support the reliable instantiation of adaptive behavior (Mount and Monje, 2017). While the role of cortical myelination in regulating plasticity during sensitive periods of neurodevelopment has been demonstrated in animal models (Takesian and Hensch, 2013), there are few in vivo studies characterizing the maturational timing of cortical myelination during sensitive periods of human development.
The cerebral cortex follows an overarching organizational axis that is anchored on one end by sensorimotor areas underlying perception and movement, and at the opposing end by heteromodal and paralimbic association areas involved with cognitive control and socioemotional processing (Huntenburg et al., 2018). The archetypal sensorimotor-association (S-A) axis of cortical organization was derived in a recent study by computing a composite ranking of cortical maps that capture feature variation in cortical microstructure, functional connectivity, metabolism, gene expression, and evolutionary expansion of the cortical mantle (Sydnor et al., 2021). Further, the spatiotemporal pattern of cortical development may also progress along the S-A axis, such that the early maturation of sensorimotor systems scaffolds the development of association areas that integrate more complex cognitive processes. Both myelin-sensitive imaging and postmortem histology have provided evidence for continuous remodeling of intracortical myelin content into adulthood (R.A. Hill et al., 2018; Hughes et al., 2018), but it remains unclear whether cortical areas exhibit graded differences from childhood through early adulthood in the developmental timing of myelination according to their position along the S-A axis.
Advances in noninvasive magnetic resonance imaging (MRI) methods have facilitated in vivo investigations of myelin content and tissue changes across the human lifespan (Edwards et al., 2018). The ratio of T1-weighted and T2-weighted MRI images (T1w/T2w) provides an indirect estimate of cortical myelin that correlates with histologic myelin stains and quantitative MRI (e.g., magnetization transfer, R1, R2) indices of cortical myelin content (Glasser and Van Essen, 2011; Glasser et al., 2014, 2021; Shams et al., 2019). T1w/T2w mapping has been increasingly used in recent years because it is based on widely available MRI sequences that are routinely collected in brain imaging studies at high resolution, at clinical field strengths, and with short acquisition times. However, prior studies evaluating age-related differences in T1w/T2w (Grydeland et al., 2019; Norbom et al., 2020) have not accounted for residual radio frequency transmit field (B1+) biases, which impact T1w/T2w estimates (as it does quantitative relaxometry measures) in an age-dependent manner through their relationship to body size (Glasser et al., 2021; MacLennan et al., 2022). Because studies that neglected B1+ bias field effects may be prone to spurious results, it is critical to characterize age-related differences in T1w/T2w after mitigating these confounds.
Here, we apply high-resolution cortical T1w/T2w mapping methods from the Human Connectome Project in Development (HCP-D; N = 628, ages 8–21) to characterize the maturational timing of cortical myelination from childhood through early adulthood. We test whether cortical areas exhibit graded differences in the rate of T1w/T2w development according to their position along the S-A axis. We use a novel approach for mitigating B1+ confound effects in T1w/T2w mapping that enables us to detect unbiased age-related differences in T1w/T2w (Glasser et al., 2021). We characterize properties of age-related differences in T1w/T2w across the cerebral cortex, including estimates of the variance in T1w/T2w explained by age, rate of change, age of peak growth, and degree of nonlinearity in T1w/T2w development from 8 to 21 years old. These data provide evidence for graded variation in age-dependent T1w/T2w differences along the S-A axis that may reflect variation in the rate of cortical myelination during adolescence and undergird the development of complex information processing.
Materials and Methods
Participants
Neuroimaging datasets included 628 typically developing participants aged 8–21 years (53.5% female; Fig. 1A) who were part of the HCP-D. The HCP-D is a large cross-sectional and longitudinal study aiming to characterize brain connectivity development in a sample approximating the demographics of youth in the United States with respect to race, ethnicity, and socioeconomic status (Harms et al., 2018; Somerville et al., 2018; Elam et al., 2021). Participants in this sample were recruited across four sites: Harvard University, University of California-Los Angeles, University of Minnesota, and Washington University in St. Louis. Exclusion criteria for recruitment included (1) premature birth (<37 weeks gestation); (2) serious neurologic condition (e.g., stroke, cerebral palsy); (3) serious endocrine condition (e.g., precocious puberty, untreated growth hormone deficiency); (4) long-term use of immunosuppressants or steroids; (5) history of serious head injury; (6) hospitalization >2 d for certain physical or psychiatric conditions or substance use; (7) treatment >12 months for psychiatric conditions; (8) claustrophobia; or (9) pregnancy or other contraindications for MRI. Participants provided written informed consent and assent, and parents of participants under 18 years provided written informed consent for their child's participation. All procedures were approved by a central Institutional Review Board administered at Washington University in St. Louis (IRB #201603135).
This study focuses on cross-sectional data drawn from the publicly available HCP-D “Release 2.0” dataset in the National Institute of Mental Health Data Archive (NDA; N = 652). From this dataset, we excluded 21 participants under the age of 8 because of insufficient data available for 5–7 year olds at the time of this study. We excluded two additional participants because of face-masking errors during image reconstruction, and we excluded one additional participant with poor T1w/T2w data quality indicated by visible artifacts and substantial outliers in T1w/T2w maps, yielding our final sample of 628 youth.
Myelin measurement: caveats and assumptions
MRI signal contrasts in cortical gray matter reflect a composite of tissue properties (Fukunaga et al., 2010; Stüber et al., 2014). Changes in the MRI signal may reflect changes in iron, myelin, lipids, water, and local cell density (Fukunaga et al., 2010; Stüber et al., 2014; Edwards et al., 2018). Several noninvasive MRI methods have been developed in recent years to estimate cortical myelin content, including quantitative T1, quantitative T2*, quantitative magnetization transfer, T2-based myelin water fraction methods (Carey et al., 2018; Edwards et al., 2018), and T1w/T2w ratio image intensities. These approaches are correlated with one another and with histologic myelin stains in the cerebral cortex (Glasser and Van Essen, 2011; Glasser et al., 2014, 2021), but the relative sensitivity to the various sources of MRI tissue contrast differ across these methods. Importantly, differences in these measures over the course of development may reflect changes in multiple tissue properties (e.g., not just myelin) and the properties that are changing may differ from those that primarily drive the overall image contrast (e.g., the known and well validated differences across the cortex between heavily myelinated and lightly myelinated cortical areas). Thus, age-dependent differences in cortical myelin measures should be interpreted cautiously as microstructural development likely reflecting myelination, while acknowledging the complexity of cortical tissue changes during development and the limitations of indirect proxy measures.
Image acquisition
High-resolution T1w MRI images were acquired on a 3T Siemens Prisma with a 32 channel head coil using a 3D multiecho MPRAGE sequence [Mugler and Brookeman, 1990; van der Kouwe et al., 2008; 0.8-mm isotropic voxels, TR/TI = 2500/1000 ms, TE = 1.8/3.6/5.4/7.2 ms, flip angle = 8°, in-plane (iPAT) acceleration factor of 2, TA = 8:22, up to 30 reacquired TRs]. Structural T2w images were acquired at 0.8 mm isotropic using the variable-flip-angle turbo-spin-echo 3D SPACE sequence (Mugler et al., 2000; TR/TE = 3200/564 ms; same in-plane acceleration, TA = 6:35, up to 25 reacquired TRs).
Both “PreScan Normalized” (Siemens' approach for removing the receive coil intensity profile) and non-normalized reconstructions were generated at the scanner. The former were used for image quality review at the scanner, while the latter were used as the inputs for subsequent processing (consistent with the use of non-normalized reconstructions as the input for the processing of the HCP Young Adult acquisitions; Van Essen et al., 2013). Both versions (of the T1w) image were used for estimating the B1– receive field to correct for effects of subject motion between the T1w and T2w images. Only the first two echoes of the T1w image were used in processing because of artifacts in the later echoes that affected surface reconstructions and T1w/T2w maps (Elam et al., 2021).
Volumetric navigators (vNavs) were embedded in the T1w and T2w sequences for prospective motion correction and for selective reacquisition of the lines in k-space that were heavily corrupted by subject motion (Tisdall et al., 2012). Real-time motion correction can substantially reduce bias in brain morphometry analyses, where motion might induce measurable morphometric differences (Reuter et al., 2015; Tisdall et al., 2016). If the T1w or T2w structural scans were nonetheless deemed to be of poor quality at the time of acquisition, they were reacquired (typically immediately in the same imaging session, although sometimes in a different session). Only the single pair of T1w and T2w scans rated the highest in quality from a single session were used in subsequent processing. For further neuroimaging protocol description, see Harms et al. (2018). Additionally, 2-mm isotropic gradient echo (GRE) and spin echo (SE) images were acquired and used for computing the pseudo-transmit field described below (Glasser et al., 2021).
Image processing
The structural MRI data were analyzed using the HCP Pipelines (Glasser et al., 2013) version 4.0.0, instantiated into the QuNex container environment (qunex.yale.edu), and the data were released as a part of the Lifespan HCP Release 2.0 in the NDA (https://nda.nih.gov/). Briefly, the T1w and T2w volumes were processed through the PreFreeSurfer pipeline, which included gradient nonlinearity distortion correction, initial brain-extraction, and rigid registration into an anterior/posterior-commissure aligned “native” space, registration of the T2w volume to the T1w volume using boundary-based registration (BBR; Greve and Fischl, 2009), correction for the receiver coil bias field based on the smoothed square root of the product of the T1w and T2w images, and registration of the structural images to MNI space. Next, the FreeSurfer pipeline (v6.0.0; Dale et al., 1999; Fischl, 2012) optimized for use with high-spatial resolution (>1 mm isotropic) structural images was used for computing the “white” and “pial” surfaces, including use of the T2w volume to optimize the pial surface placement. Lastly, the PostFreeSurfer pipeline produced cortical surface models in GIFTI format and surface-related data in CIFTI format, and each subject's cortical surface was then registered to a common 32k_FS_LR mesh using “MSMAll” areal-feature-based cortical surface registration, which is a multimodal registration constrained by cortical T1w/T2w maps and resting-state network maps (Robinson et al., 2014, 2018; Glasser et al., 2016b).
Following cortical surface reconstruction, a single experienced individual performed a “SurfaceQC” review of the white and gray matter surface placement, informed by the T1w/T2w ratio maps (Glasser and Van Essen, 2011; Elam et al., 2021). Participants with more than minor (focal) issues were flagged for possible future editing and excluded from the cohort analyzed for the current study. This “SurfaceQC” review of the HCP-D data revealed some degradation of the accuracy of surface placement relative to expectations established by the HCP Young Adult project, which we subsequently traced to artifacts in the longer echoes. Therefore, to reduce the prevalence of surface segmentation errors in this developmental sample, we used the mean of the shortest two echoes (i.e., excluded the longest two of four echoes) as the T1w input to the HCP Pipelines (Elam et al., 2021).
T1w/T2w processing
T1w/T2w maps were generated as a cortical ribbon volume and a surface map using the methods described previously (Glasser and Van Essen, 2011; Glasser et al., 2013), together with more recent improvements (Glasser et al., 2014, 2021). Briefly, the T1w/T2w ratio in the voxels between the white and pial surfaces is mapped to the surface mesh in a way that emphasizes the middle layers and underemphasizes voxels near pial and gray/white surface to reduce partial volume effects (Glasser et al., 2013). Division of the T1w image by the T2w image mathematically cancels the signal intensity bias related to the sensitivity profile of the radio frequency receiver coils, which is the same in both images in the absence of subject head motion. Taking the ratio will also increase the contrast related to myelin content since both input images show myelin-related contrast (Glasser and Van Essen, 2011), inverted in the T2w image relative to the T1w image. Individual T1w/T2w maps were parcellated according to the HCP-multimodal atlas (Glasser et al., 2016a) using wb_command -cifti-parcellate in Connectome Workbench v1.4.2 (Marcus et al., 2011). This computes the average value of all vertices within a parcel and is thus a neuroanatomically constrained form of smoothing. T1w/T2w units are arbitrary, representing relative estimates of intracortical myelin content that are comparable within a consistently acquired study. Further, since the scaling of T1w and T2w images depends on the scanner and acquisition parameters, T1w/T2w units should not be directly compared across studies without appropriate harmonization.
B1+ transmit field correction of T1w/T2w myelin maps
T1w/T2w maps were initially developed for neuroanatomical analyses such as constraining the anatomic boundaries of cortical areas within individuals (Glasser et al., 2016a) and residual biases from the B1+ radio frequency transmit field were mitigated using the MyelinMap_BC approach (Glasser et al., 2013); however, this approach also removes genuine low spatial frequency individual differences in the T1w/T2w ratio. Given the increasing interest in using T1w/T2w maps to perform statistical comparisons across individuals and groups with other variables of interest, such as age (Kwon et al., 2018; Grydeland et al., 2019; Norbom et al., 2020), a new approach to B1+ transmit field bias correction was developed (Glasser et al., 2021). B1+ biases are influenced by the loading of the body coil (influenced by both the participant's head and body). Thus, variables such as head size, body size, and body mass index (BMI), which are often correlated with age and other variables of interest, can modulate B1+ bias and may lead to potentially spurious results in cross-subject statistical analyses of T1w/T2w maps.
In this study, we applied a novel, empirically validated “pseudo-transmit field” correction to mitigate B1+ bias in individual T1w/T2w maps, thereby reducing potentially spurious age-related differences in T1w/T2w values (Glasser et al., 2021). The B1+ correction approach in this study relies on computing a pseudo-transmit field based on the average across phase encoding directions of GRE images divided by SE images (mimicking the GRE T1w scan divided by SE T2w scan). First, a reference T1w/T2w map was generated at the group level by finding the scaling between the group average pseudo-transmit field and group average T1w/T2w map that minimizes the correlated left-right differences between the two maps (i.e., the clearly spurious left-right asymmetries). This reference group map was used to correct the individual maps. For the individual correction, the pseudo-transmit map was scaled to minimize the correlated differences between the individual's T1w/T2w map and the reference T1w/T2w map and the pseudo-transmit map (which includes all differences, not simply left-right ones, and is more robust at the individual level). Before estimating this correction, any residual B1– effects because of subject head motion between the T1w and T2w images were also removed using the scanner-computed B1– receive field. The pseudo-transmit field requires regularization by thresholding regions of T2*-related signal loss combined with spatial smoothing (with compensation for intensity changes induced by smoothing); it is then scaled to equal 1 at the value where the GRE/SE ratio corresponds to the flip angle prescribed by the scanner, a reference value that is determined at the group level. For more methodological details and validation of the correction, see Glasser et al. (2021), which demonstrates that this correction approach eliminates spurious relationships between the T1w/T2w ratio and B1+ bias modulated by body size.
Statistical analysis
Bayesian generalized additive models were fitted using the brms (Bürkner, 2018) interface to the Stan modeling language (Stan Development Team, 2021). Models were fitted separately to each parcellated cortical area to characterize different properties of age-related changes in T1w/T2w across the cerebral cortex. This approach samples from the posterior smooth function of T1w/T2w conditional on age and covariates, representing age-related changes in T1w/T2w, and the posterior derivative of the smooth function on age, representing the rate of change across the age range. Importantly, Bayesian thin-plate spline models allow us to estimate both linear and nonlinear age-related changes in T1w/T2w without specifying a linear or polynomial function a priori (Fahrmeir and Kneib, 2011; Hastie and Tibshirani, 2017; Wood, 2017). All statistical analyses were conducted in R version 3.5.1 (R Core Team, 2018).
We evaluated age-related change in T1w/T2w for each cortical parcel while controlling for participant sex, scanner, and several covariates related to the B1+ transmit field correction including the scanner transmit voltage, the mean of the pseudotransmit map, four regularization parameters [(1) T2* dropout threshold; (2) smoothing FWHM; (3) correction factor for smoothing's effect on the pseudotransmit field's intensities; and (4) the slope parameter of the correction], and a corrected T1w/T2w lateral ventricular CSF regressor (after excluding partial volume voxels and CSF flow effects) that are described in detail previously (Glasser et al., 2021). The brms syntax for the full Bayesian model for each cortical parcel was expressed as follows:
brms(T1w/T2w ∼ s(Age, k = −1, bs = 'tp') + Sex + Scanner + Transmit_Voltage + Mean_Pseudo_Transmit_Map + T2_Dropout_Threshold + Smoothing_FWHM_mm + Smoothing_Correction_Parameter + Correction_Slope + Corrected_CSF_Regressor)
Thin plate regression splines were used for the smoothing basis, and k was left at the default value of −1 which entails setting the basis dimension to 10. We retained the default priors set by brms; for regression coefficients this is a flat prior over the reals; for the intercept this is a Student's t distribution with df = 3, µ set at the mean of the outcome, and σ = 2.5; for the SD of the splines and SD of the error term this is a Student's t distribution with df = 3, µ = 0, and σ = 2.5. We ran four chains with 4500 total iterations (2000 warm-up) for each chain, yielding 10,000 posterior draws in total. The posterior smooth function for each cortical parcel was used to estimate other properties of age-related change in T1w/T2w, as described below.
Measurements of age-related change in T1w/T2w
Partial R2 values for age splines in each regional model were calculated by taking the difference between full models including an age spline and reduced models without age. These R2 values reflect the amount of variance in cortical T1w/T2w explained by age.
Annualized rate of change (AROC) in T1w/T2w
The annualized estimated rate of change in T1w/T2w was calculated for each cortical parcel by taking the difference between posterior samples for T1w/T2w at age 21 and age 8, and then dividing that difference by the number of years in the age range of the sample. This yielded a posterior distribution for this difference which we use to compute a point estimate using the median and 95% credible interval. The estimated rate of change in T1w/T2w was based on cross-sectional differences and does not reflect direct longitudinal estimates.
Linearity of age-related changes in T1w/T2w
Given prior work showing variation in the linearity or shape of neurodevelopmental changes during adolescence (Shaw et al., 2008), we used two complementary approaches to evaluate the linearity of age-related differences in T1w/T2w for each cortical parcel. Both approaches were dependent on the Bayesian spline models described above, which estimate both linear and nonlinear age-related differences in T1w/T2w without specifying a linear or polynomial function a priori (Wood, 2017). Therefore, to estimate the linearity of age curves regardless of functional form, we defined a measure of linearity as the mean absolute posterior second derivative of the smooth function of T1w/T2w on age.
First, the mean absolute posterior second derivative of the smooth function of T1w/T2w on age was calculated as a continuous measure of linearity for each cortical parcel using the curvish package in R (Flournoy, 2021). This measure reflects the rate of change in the slope of T1w/T2w age curves (e.g., how the rate of T1w/T2w change varies from 8 to 21 years).
Second, we conducted a separate, more conservative test evaluating whether a (nonlinear) spline model or a linear model of age-related changes in T1w/T2w provides greater out-of-sample predictive performance based on the leave-one-out-information criterion (LOOIC; Vehtari et al., 2017). Specifically, the smooth function from a generalized additive model was identified as nonlinear when the LOOIC difference between a Bayesian linear and spline model was at least 2.0 times the approximate standard error of its estimate (Bürkner, 2018). This is approximately analogous to a two-sided Z-test. The LOOIC indexes the out-of-sample predictive performance of a model. Therefore, cortical parcels for which the generalized additive model fits better are those in which we have credible evidence that this more complex, nonlinear form has greater out-of-sample predictive performance.
Age of peak and slowing T1w/T2w growth
The posterior age of peak growth (i.e., steepest slope) in T1w/T2w was estimated using the curvish package in R (Flournoy, 2021). Using the posterior of the smooth function and the method of finite differences at 100 equally spaced points, we evaluate the first derivative. Then, for each posterior sample we find the age at which the derivative is maximized, yielding a posterior density reflecting the probability of the slope being maximized at any given age.
To identify windows in development where the rate of age-related increases in T1w/T2w credibly slowed down, we computed the difference between the posterior derivative of the smooth at each evaluation point and the posterior derivative at the age of maximum median derivative. Regions of this 95% credible interval that did not include zero were used to mark regions of the smooth that were credibly less steep than the slope at the point of greatest maturation.
Identifying cortical areas with similar rate and shape of age-dependent T1w/T2w changes using clustering analysis
Based on prior literature characterizing variation in the maturational timing of cortical neurodevelopment across brain areas (Shaw et al., 2008; Somerville, 2016; Sydnor et al., 2021), we aimed to identify data-driven clusters of cortical parcels with similar T1w/T2w smooth functions on age. Further, rather than clustering parcels based on a single scalar property, we were interested in clustering cortical parcels based on latent dimensions of age-dependent T1w/T2w changes that capture subtle variation in the rate and shape of cortical development. To this end, we applied a functional latent mixture model using the funHDDC package in R (Bouveyron and Jacques, 2011), which consistently outperforms alternative clustering algorithms applied on spline coefficients (Bouveyron and Jacques, 2011; Schmutz et al., 2020).
First, for Bayesian models in each cortical parcel, we subtracted the posterior T1w/T2w estimates at age 8 from predicted T1w/T2w values at 500 age bins distributed equally between 8 and 21 years. This effectively removed intercept differences (i.e., starting T1w/T2w level) across cortical parcels, which were not of interest in this analysis. Second, we converted a brain region by age matrix into a functional data object for input into the funHDDC algorithm. Third, we ran funHDDC on the functional T1w/T2w trajectories, allowing the optimal number of clusters to vary freely between 2 and 7. Fourth and finally, we used the funHDDC::slopeHeuristic function to identify the most parsimonious model solution based on penalized log-likelihood (Schmutz and Bouveyron, 2021).
Sensitivity analyses
Cortical thickness
Prior studies have shown that T1w/T2w and cortical thickness measures are inversely correlated across the cerebral cortex (Glasser and Van Essen, 2011; Glasser et al., 2014; Natu et al., 2019). To demonstrate that our results were specific to age-related differences in T1w/T2w and are not attributable to concurrent changes in cortical thickness, we fit Bayesian spline models for each cortical parcel while controlling for estimates of cortical thickness (i.e., “*_V1_MR.thickness.32k_fs_LR.dscalar.nii”). We then calculated the spatial correlation between parcellated cortical maps of the AROC in T1w/T2w with and without controlling for cortical thickness.
Head motion
Head motion during structural MRI acquisition can impact image quality, inducing spurious artifacts in image intensity and morphometric measurements (Tisdall et al., 2016; Savalia et al., 2017). Age-related differences in motion (younger individuals commonly move more than older individuals) are a common concern when evaluating age-related changes in brain imaging measures (Satterthwaite et al., 2013; Reuter et al., 2015; Baum et al., 2018). Thus, we undertook sensitivity analyses in which participants with the highest motion estimates were excluded to test whether the primary results were evident in this lower-motion sample. To do so, we estimated the AROC in T1w/T2w after excluding 114 participants who exceeded the number of allowed TR reacquisitions for either the T1w (up to 30 allowed) or T2w (up to 25 allowed) acquisition, and thus had a considerable degree of ongoing head motion that was not fully corrected. For these participants, the number of k-space lines that still exceeded the motion threshold after the allowed number of reacquisitions is unknown, and thus the quality of these images is degraded relative to what could have been achieved if there was no limit placed on the number of reacquisitions (Tisdall et al., 2012). We then calculated the spatial correlation between parcellated cortical maps of the AROC in T1w/T2w with and without including participants with excess head motion.
Linking T1w/T2w development with the S-A axis of cortical organization
To evaluate whether cortical parcels exhibited differences in age-related changes in T1w/T2w according to their position along the S-A axis, we calculated the Spearman correlation between parcellated maps of cortical T1w/T2w development and S-A axis rankings. We used a spatial permutation test (described in detail below) to assess the significance of spatial correlations between each measure of age-related change in T1w/T2w described above (partial R2, rate of change, age of peak growth, degree of nonlinearity) and the S-A axis of cortical organization.
Briefly, S-A rankings were derived from cortical maps of ten fundamental brain features exhibiting systematic variation between lower-order sensorimotor areas and higher-order association areas (Sydnor et al., 2021). These cortical maps measured the anatomic hierarchy from T1w/T2w mapping (Glasser and Van Essen, 2011), the functional hierarchy quantified by the principal gradient of functional connectivity (Margulies et al., 2016), the evolutionary hierarchy quantified by macaque-to-human cortical expansion (J. Hill et al., 2010), allometric scaling quantified as the relative extent of areal scaling with overall brain size (Reardon et al., 2018), aerobic glycolysis quantified from positron emission tomography measures of oxygen consumption and glucose use (Vaishnavi et al., 2010), cerebral blood flow quantified by arterial spin labeling (Satterthwaite et al., 2014), gene expression quantified by the first principal component of brain-expressed genes (Burt et al., 2018), the first principal component of NeuroSynth meta-analytic decodings (Yarkoni et al., 2011), externopyramidization quantified as the ratio of supragranular pyramidal neuron soma size to infragranular pyramidal neuron soma size (Paquola et al., 2020b), and cortical thickness quantified from structural MRI (Human Connectome Project S1200 data). See Sydnor et al. (2021) for further details.
Spatial permutation test
The statistical significance of spatial correlations (Pspin) between parcellated cortical maps of age-related change in T1w/T2w and the archetypal S-A axis were assessed using a conservative parcel-based spatial permutation test that preserves the spatial covariance structure of cortical maps (Alexander-Bloch et al., 2018), as implemented by Váša et al. (2018). This permutation procedure generates a distribution of 10,000 spatial null maps randomly rotated on the spherical cortical surface to retain spatial contiguity and hemispheric symmetry of the original cortical maps. For each parcellated map showing age-related change in T1w/T2w (i.e., R2, rate of change, age of peak growth, nonlinearity, latent subgroups), we generated 10,000 randomly rotated “null” maps on the cortical surface. We then estimated the null distribution of Spearman correlation coefficients between the archetypal S-A axis map and the randomly rotated maps of T1w/T2w development. Permutation-based p-values (Pspin) were calculated based on the number of times the empirical correlation coefficient was higher than the permuted correlation coefficient.
False-positive error rate correction
When evaluating statistical significance of spatial correlations between cortical maps of age-related change in T1w/T2w and the archetypal S-A axis, the spatial permutation (“spin”) test provides family-wise control of multiple comparisons (Type I error) by accounting for the spatial and contralateral dependence of areal measurements in cortical maps (Alexander-Bloch et al., 2018). Further, we used the Holm correction (Holm, 1979), which is more powerful than Bonferroni and valid under arbitrary assumptions, to adjust the p-values across the 5 tests that are used to test the hypothesis that T1w/T2w development is linked with the S-A axis of cortical organization. The same correction was applied across all 15 spatial correlation tests between pairs of parcellated maps showing different properties of T1w/T2w development and the S-A axis (α = 0.05; Table 1). For descriptive purposes, we dichotomize parcels into those with or without a more complex functional form, and those with or without evident windows where the rate of age-related T1w/T2w growth credibly slows down, comprising 360 tests for each set. We do not attempt to control the family-wise false positive rate for the linearity decisions, or probability of making a sign error (Gelman and Carlin, 2014) when identifying windows where the rate of age-related T1w/T2w growth credibly slows down, as we report these as hypothesis-generating descriptions for which a higher Type I error rate (and lower Type II error rate) is preferable.
Data access and availability
Minimally preprocessed structural neuroimaging data from the HCP-D is available for download through the NIMH Data Archive (https://nda.nih.gov/). Parcellated cortical maps presented in this study are available for download through the BALSA database (https://balsa.wustl.edu/study/P2DmK).
Results
We analyzed high-resolution cortical T1w/T2w maps from 628 participants (ages 8–21 years old; Fig. 1A,B) to characterize the development of cortical T1w/T2w during youth after correcting for transmit field biases (Glasser et al., 2021). The posterior smooth function of T1w/T2w on age and the posterior derivative of the age function were estimated by fitting Bayesian thin-plate spline models across 180 cortical parcels in each hemisphere, thereby obtaining measures of the overall effect size, rate, and nonlinearity of age-related change in T1w/T2w across the neocortical sheet. For example, we illustrate age-related differences in T1w/T2w estimated for the left primary motor cortex using the posterior smooth function (Fig. 1C) and the posterior derivative (Fig. 1D) of T1w/T2w on age.
Charting T1w/T2w development during youth. A, Age/sex histogram of the current sample of 628 youth who completed structural neuroimaging as part of the HCP-D. B, T1w/T2w maps were parcellated using the HCP multimodal atlas (Glasser et al., 2016a) and averaged across participants. T1w/T2w units are arbitrary, representing relative estimates of intracortical myelin content that are comparable within a consistently acquired study. T1w/T2w units of 1.4 and 1.9 correspond to the second and 98th percentiles in this dataset. Data from the left primary motor cortex (highlighted with black border) are shown in panels C, D. C, We fit Bayesian generalized additive models with thin-plate splines to estimate different properties of age-related change in cortical T1w/T2w. Specifically, we estimated the posterior smooth function of T1w/T2w on age for each of the 360 cortical parcels (data shown for the left primary motor cortex, highlighted in panel B). The shaded area represents the 95% credible interval of the posterior smooth function. Navy blue segments of the posterior smooth function indicate the slope is credibly as steep as the maximum slope, while yellow segments of the posterior smooth function indicate the slope is credibly less steep than the maximum slope (see panel D). D, Bayesian generalized additive models were fitted to estimate the posterior derivative of the smooth function of T1w/T2w on age to characterize the rate of change, with higher values indicating a steeper slope of change per unit age and 0 representing a flat slope (i.e., no age-related change). The shaded area represents the 95% credible interval of posterior derivative estimates. The vertical line marks the age of the maximum median derivative (steepest slope). To identify windows where the rate of age-related T1w/T2w growth credibly slows down, we computed the difference between the posterior derivative of the smooth and the posterior derivative at the age of the maximum median derivative; regions of this 95% credible interval (data not shown) that did not include zero were used to mark regions of the smooth that were credibly less steep than the slope at the point of greatest maturation. Data are shown for an exemplar cortical parcel to illustrate our analytical approach.
For a descriptive visualization of cortical T1w/T2w development, we averaged data into three age groups (8–10 years approximating childhood, 14–16 years approximating mid-adolescence, and 19–21 years approximating emerging adulthood) to illustrate the spatial patterning of age-related differences in T1w/T2w from childhood to adulthood. The pattern of cortical T1w/T2w differences during adolescent development reflects the S-A axis reported by Sydnor et al. (2021; Fig. 2, inset right). The S-A hierarchy was evident in 8–10 year olds (n = 145; Fig. 2, left), with high T1w/T2w in sensory cortex and low T1w/T2w in association cortex. This hierarchy was present across all age groups, suggesting that the S-A hierarchy in cortical T1w/T2w is reinforced from childhood through adulthood. Relative to younger participants, 14–16 year olds (n = 154; Fig. 2, middle) had higher T1w/T2w across sensorimotor and intermediate association areas of the cerebral cortex. In 19–21 years olds, we observed a spatial pattern of cortical T1w/T2w extending from primary sensory to multimodal association cortex (n = 120; Fig. 2, right).
Waves of cortical T1w/T2w maturation from childhood through early adulthood. Parcellated T1w/T2w maps were averaged across participants within three age groups to illustrate the spatial patterning of age-related increases in T1w/T2w during youth. T1w/T2w units of 1.4 and 1.9 correspond to the second and 98th percentiles in this dataset. Age groups include 8–10 year olds (left; n = 145), 14–16 year olds (middle; n = 154), and 19–21 year olds (right; n = 120). Age-related increases in cortical T1w/T2w were observed across the cortical sheet, reinforcing a S-A hierarchy in cortical microstructure. The archetypal S-A axis (inset right) was adapted with permission from Sydnor et al. (2021).
Next, we evaluated the topography of the effect size and rate of age-related change in T1w/T2w from childhood through adulthood. In a multiple linear regression, all cortical parcels showed significant age-related differences in T1w/T2w after family-wise error correction (p < 0.05) with the exception of the bilateral subgenual cingulate (Brodmann area 25), which suffers from low SNR because of proximity of the ventricles and medial wall. Moreover, we found systematic variation in the amount of variance in T1w/T2w explained by age. Specifically, age explained a greater proportion of variance in T1w/T2w in sensorimotor areas such as right area V2 (R2 = 0.25; Fig. 3B) compared with association zones where age explained little variance in T1w/T2w, such as the left anterior cingulate cortex (area p32; R2 = 0.028; Fig. 3C). Figure 3D shows the spatial correlation between the variance in T1w/T2w explained by age (R2) and the cortical parcel's position along the archetypal S-A axis of cortical organization (r = −0.65, Pspin < 0.001).
Effect size estimates for regional Bayesian models of T1w/T2w development. A, Partial R2 values for age splines were estimated for each cortical parcel, representing the proportion of variance in T1w/T2w differences explained by age. B, Age-related increases in T1w/T2w were relatively strong in sensorimotor areas such as the right V2 (outlined and labeled in panel A), where age explained 25% of variance in T1w/T2w. C, By contrast, age-related increases in T1w/T2w were relatively weak in heteromodal and paralimbic association areas such as the left anterior cingulate cortex (area p32; outlined and labeled in panel A), where age explained 2.8% of variance in T1w/T2w. D, The variance in T1w/T2w explained by age (R2) was correlated with the cortical area's position along the archetypal S-A axis of cortical organization. In panel D, the blue data point corresponds to the exemplar sensorimotor parcel (V2) and the yellow data point corresponds to the exemplar association parcel (ACC) highlighted in panels A–C. ACC = anterior cingulate cortex; S-A = sensorimotor-association. Pspin is the permutation-based p-value calculated from a conservative parcel-based spatial permutation (“spin”) test. Data are shown for exemplar sensorimotor and association parcels to illustrate differences in T1w/T2w development. Models were estimated independently for all 180 cortical parcels in each hemisphere (360 in total).
We also observed spatial variation in the AROC in T1w/T2w estimated for each cortical parcel (Fig. 4A). Sensorimotor areas had relatively high AROC in T1w/T2w. Figure 4B shows the posterior smooth for the right motor cortex (area 4), an exemplar sensorimotor area with relatively high AROC in T1w/T2w (AROC = 0.025). Figure 4C shows the posterior smooth for the left anterior cingulate cortex (area p32), an exemplar association area with relatively low AROC in T1w/T2w (AROC = 0.005). From childhood through early adulthood, spatial variation in the AROC in T1w/T2w reflected S-A hierarchy (r = −0.65, Pspin < 0.001; Fig. 4D), with heavily myelinated sensorimotor areas showing rapid T1w/T2w increase relative to heteromodal and paralimbic association areas, which had slower age-related increases in T1w/T2w.
Annualized rate of change in T1w/T2w during youth. A, Annualized rate of change (AROC) in T1w/T2w was estimated for each cortical parcel. B, The AROC was relatively high in sensorimotor areas such as the right primary motor cortex (area 4; highlighted in panel A, right). C, In contrast, the AROC in T1w/T2w was relatively low in prefrontal and paralimbic association areas such as the left medial prefrontal cortex (area 9m; highlighted in panel A, left). D, The AROC in T1w/T2w was correlated with the cortical parcel's position along the archetypal S-A axis of cortical organization. The blue data point corresponds to the exemplar sensorimotor parcel (primary motor cortex) and the yellow data point corresponds to the exemplar association parcel (medial PFC) highlighted in panels A–C. Shaded areas in the plots of panels B, C represent the 95% credible interval of the posterior smooth function on age. Segments of the age curve where the slope is credibly less steep than the maximum slope (indicating a credibly reduced rate of change) are highlighted in yellow. PFC = prefrontal cortex; S-A = sensorimotor-association. Pspin is the permutation-based p-value calculated from a conservative parcel-based spatial permutation (“spin”) test. Data are shown for exemplar sensorimotor and association parcels to illustrate differences in T1w/T2w development. Models were estimated independently for all 180 cortical parcels in each hemisphere (360 in total).
Sensitivity analyses revealed high correspondence between the original parcellated cortical map of AROC in T1w/T2w and the parcellated map when controlling for cortical thickness in each cortical area (r = 0.96, Pspin < 0.001). We also observed a highly consistent spatial pattern of T1w/T2w development when excluding participants with excess head motion (n = 451; r = 0.97, Pspin < 0.001).
Next, we evaluated whether cortical parcels varied systematically in the linearity of age-related increases in T1w/T2w, as nonlinear T1w/T2w smooth functions might reflect diminishing age-related changes or a developmental inflection point. Figure 5A shows the linearity of age-related differences in T1w/T2w using the mean absolute posterior second derivative, which captures change in the slope over time. Figure 5B shows the posterior smooth function of T1w/T2w on age for the left V2, an exemplar sensorimotor area with relatively nonlinear age-related increases in T1w/T2w (mean absolute second derivative = 0.004). Figure 5C shows the posterior smooth function of T1w/T2w on age for the right anterior ventral insula, an exemplar association area with relatively linear age-related increases in T1w/T2w (mean absolute second derivative = 0.002). Spatial variation in the linearity of age-related increases in T1w/T2w was significantly correlated with the cortical parcel's position along the S-A axis (r = −0.45, Pspin < 0.001; Fig. 5D), with sensorimotor areas showing relatively nonlinear T1w/T2w increase relative to heteromodal and paralimbic association areas, which had relatively linear T1w/T2w increase with age. The LOOIC difference between linear and nonlinear spline models indicated that 63 of 360 cortical parcels (17.5%) had better out-of-sample predictive performance when including the spline. We found that 116 of 360 (32.2%) cortical parcels had segments of the age curve that were credibly less steep than the maximum slope, indicating a credible decrease in the slope of T1w/T2w maturation over time. This deceleration of age-related increases in T1w/T2w occurred primarily between 18 and 21 years old.
Nonlinearity of age-related changes in T1w/T2w during youth. A, The nonlinearity of age-related changes in T1w/T2w was estimated for each cortical parcel using the mean absolute posterior second derivative, where higher values indicate more nonlinear growth. B, Sensorimotor areas such as the left V2 (highlighted in panel A, left) had relatively nonlinear T1w/T2w growth, with the rate of change decreasing credibly by 17.5 years old (indicated by yellow segment). C, Cortical parcels in heteromodal and paralimbic association areas such as the right anterior ventral insula (highlighted in panel A, right) had relatively linear growth in T1w/T2w (constant slope). D, The nonlinearity of age-related changes in T1w/T2w was correlated with the cortical parcel's position along the archetypal S-A axis of cortical organization (Sydnor et al., 2021). Sensorimotor areas exhibited significantly higher nonlinearity in age-related increases in T1w/T2w compared with areas of association cortex. In panel D, the blue data point corresponds to the exemplar sensorimotor parcel (V2) and the yellow data point corresponds to the exemplar association parcel (anterior insula) highlighted in panels A–C. Shaded areas in the plots of panels B, C represent the 95% credible interval of the posterior smooth. Segments of the posterior smooth where the slope is credibly less steep than the maximum slope (indicating a credibly reduced rate of change) are highlighted in yellow. Pspin is the permutation-based p-value calculated from a conservative parcel-based spatial permutation (“spin”) test. Data are shown for exemplar sensorimotor and association parcels to illustrate differences in T1w/T2w development. Models were estimated independently for all 180 cortical parcels in each hemisphere (360 in total).
We further characterized cortical variation in the maturational timing of peak age-related increases in T1w/T2w based on the median posterior age where the first derivative was maximized (i.e., age of steepest slope; Fig. 6A). Figure 6B shows the shape and posterior age of peak T1w/T2w growth for the left V1, an exemplar sensorimotor area with a relatively early age of peak growth in T1w/T2w (median age = 9.6 years). Figure 6C shows the shape and posterior age of peak T1w/T2w growth for the right prefrontal area 8c, an exemplar association parcel with a relatively late age of peak growth in T1w/T2w (median age = 15.0 years). Spatial variation in the age of peak T1w/T2w growth from 8 to 21 years old closely followed the S-A axis of cortical organization (r = 0.60, Pspin < 0.001; Fig. 6D), with sensorimotor areas showing relatively early peak growth in T1w/T2w compared with heteromodal and paralimbic association areas, which had relatively late peak T1w/T2w growth. In general, areas with more linear T1w/T2w growth had greater uncertainty in the age of maximum slope, as illustrated by the width of the posterior distribution (Fig. 6B,C, upper inset). This may account for why the scatter (deviation from the best-fitting line) increases with S-A ranking in Figure 6D.
Age of peak T1w/T2w growth during youth. A, The age of peak growth (i.e., steepest slope) in T1w/T2w myelin was estimated for each cortical parcel as the median posterior age where the slope is at its maximum. B, Sensorimotor areas such as left V1 (highlighted in panel A) had relatively early age of peak T1w/T2w myelin growth (median age = 9.6 years). C, Association areas such as the right prefrontal cortex (area 8c; highlighted in panel A) had relatively later age of peak T1w/T2w myelin growth (median age = 15.0 years). D, The age of peak T1w/T2w myelin growth was correlated with the cortical parcel's position along the archetypal S-A axis of cortical organization. Sensorimotor areas exhibited a significantly earlier age of peak growth compared with areas of association cortex. In panel D, the blue data point corresponds to the exemplar sensorimotor parcel (V1) and the yellow data point corresponds to the exemplar association parcel (PFC) highlighted in panels A–C. Segments of the age curve where the slope is credibly less steep than the maximum slope (indicating a credibly reduced rate of change) are highlighted in yellow. Insets at the top of panels B, C show the posterior density distribution of the age of maximum slope for exemplar areas. The shaded blue area represents the 95% credible interval of the posterior distribution; black vertical lines mark the median of the posterior distribution (values projected on the cortical surface in panel A). PFC = prefrontal cortex. Pspin is the permutation-based p-value calculated from a conservative parcel-based spatial permutation (“spin”) test. Data are shown for exemplar cortical parcels to illustrate differences in T1w/T2w development. Models were estimated independently for all 180 cortical parcels in each hemisphere (360 in total).
We assessed whether groups of cortical parcels demonstrated similar temporal aspects of age-related change in T1w/T2w, as such groups of cortical parcels might reflect the development of coordinated patterns of brain connectivity. A functional latent mixture model was used to identify groups of cortical parcels exhibiting similar age-related patterns of T1w/T2w growth during youth. The most parsimonious model identified three latent clusters or groups of cortical parcels with graded differences in the rate of age-related changes in T1w/T2w (Fig. 7A). The spatial patterning of the functional latent clusters closely recapitulated the S-A axis (Sydnor et al., 2021; F = 166.6, p < 1.0 × 10−10): sensorimotor areas at one end of the S-A axis were grouped in latent cluster 1, intermediate association areas were grouped in latent cluster 2, and heteromodal/paralimbic association areas at the opposite end of the S-A axis were grouped in latent cluster three based on variation in the developmental timing of T1w/T2w changes (Fig. 7B). Cluster 1 (Fig. 7C, magenta) primarily included sensorimotor parcels, which had the highest rate of cortical T1w/T2w growth. Cluster 2 (Fig. 7C, orange) included heteromodal association parcels showing an intermediate rate of cortical T1w/T2w growth. Cluster 3 (Fig. 7C, yellow) included paralimbic association parcels in medial prefrontal and insular cortex, which had the lowest rate of cortical T1w/T2w growth.
Graded variation across data-driven clusters of cortical T1w/T2w development. A functional latent mixture model was used to identify groups of cortical parcels exhibiting similar age-related patterns of T1w/T2w growth during youth. A, The most parsimonious model identified three latent clusters or groups of cortical parcels with graded differences in the rate of age-related changes in T1w/T2w. The inset shows the penalized log-likelihood for each model solution from k = 2 to k = 7 clusters. The best fitting three-cluster solution is highlighted in red. B, Functional latent clusters closely aligned with the S-A axis of cortical organization. Sensorimotor areas, primarily in cluster 1 (magenta), had the highest rate of cortical T1w/T2w growth. We observed graded decreases in the rate of cortical T1w/T2w growth in cluster 2 (orange), which included heteromodal association parcels in frontoparietal and temporal cortex, and further in cluster 3 (yellow), which included paralimbic association parcels in the medial prefrontal and insular cortices. C, Regional cluster assignments based on T1w/T2w development reflect the S-A axis. L-L = log-likelihood.
We observed similar spatial patterning across each measurement of age-related changes in T1w/T2w. Table 1 shows the Spearman correlation coefficients and permutation-based p-values (Pspin) for each pairwise correlation among measurements. The high correspondence between these measures suggests that different temporal properties of T1w/T2w changes are yoked during adolescence and depend on the cortical parcel's position along the S-A axis.
Correlated properties of T1w/T2w development and cortical S-A organization
Discussion
Here, we characterize the developmental timing of microstructural changes across the cerebral cortex during adolescence using recent advances in high-resolution T1w/T2w mapping applied to data from the HCP-D. This study is the first to evaluate the maturational timing of estimated cortical myelination in youth using B1+ transmit field-corrected T1w/T2w myelin maps, providing unbiased estimates of age-related differences in cortical microstructure. We demonstrate graded variation across the cortex in the timing of T1w/T2w changes during youth, with rapid age-related increases in sensorimotor cortex, and more gradual age-related increases in association cortex. This spatial pattern of microstructural brain development recapitulates the sensorimotor-to-association axis of cortical organization and plasticity during ontogeny. Moreover, our findings provide new evidence for protracted T1w/T2w changes in distributed association zones during youth, which could reflect ongoing cortical myelination underlying complex information processing and psychological functioning.
Graded variation in cortical T1w/T2w development during youth reflects the S-A axis of cortical organization
We observed graded variation in the rate and timing of age-related increases in cortical T1w/T2w from 8 to 21 years old. For each cortical parcel, Bayesian spline models provided posterior estimates of the AROC, the age of peak T1w/T2w growth, and the degree of nonlinearity in age-related changes. Each of these temporal properties of T1/T2w development followed a cortical topography reflecting S-A organization (Sydnor et al., 2021). Sensorimotor areas at one end of the S-A axis had a high rate of T1w/T2w growth, with early peak growth and relatively nonlinear age-related increases in T1w/T2w. Heteromodal and paralimbic association areas at the opposing end of the S-A axis had a slower rate of T1w/T2w changes, with late peak growth and relatively linear age-related increases in T1w/T2w. This coupling between measurements of age-related change in T1w/T2w suggest that the microstructural differentiation of cortical areas follows distinct developmental programs depending on the area's position along the S-A axis of cortical organization.
A functional latent mixture modeling procedure provided data-driven validation of the systematic variation we observed in T1w/T2w development. The most parsimonious clustering solution identified three clusters of cortical parcels based on the rate of age-dependent differences in T1w/T2w. This approach confirmed that the rate of cortical T1w/T2w growth during youth parallels the areal position along the S-A axis.
Our results provide new in vivo evidence aligned with histologic studies demonstrating topographic variation in the rate of cortical gray matter myelination during childhood and adulthood (Flechsig, 1901; Yakovlev and Lecours, 1967), underscoring a neurodevelopmental S-A axis. Graded variation in the rate of cortical myelination could reflect cascading sensitive windows for refining increasingly complex cognitive functions (Takesian and Hensch, 2013). Our findings are also consistent with comparative postmortem studies demonstrating prolonged myelination in human association cortex relative to chimpanzees (Miller et al., 2012), which could reflect an uniquely extended window of ongoing plasticity in cortical circuits underlying complex socioemotional and cognitive processing (Larsen and Luna, 2018).
Our main findings also converge with prior neuroimaging studies evaluating age-related differences in cortical myelin content using quantitative MRI methods such as magnetization transfer (MT) imaging, which demonstrate relatively late myelination in association cortices during adolescence into the third decade of life (Whitaker et al., 2016; Ziegler et al., 2019; Paquola et al., 2020a). Our findings are also consistent with a recent study evaluating age-related changes in T1w/T2w from 3–21 years old, which identified global increases in T1w/T2w across the neocortex (Norbom et al., 2020). Further, a lifespan study of 8–85 year olds characterized an early age of peak T1w/T2w growth in sensorimotor areas, and a later age of peak T1w/T2w growth in association areas (Grydeland et al., 2019). Despite their strengths, these studies did not control for B1+ transmit field bias, which impacts T1w/T2w estimates in a manner that is systematically related to age, sex, and body size (Glasser et al., 2021; MacLennan et al., 2022). Our study extends prior work by rigorously characterizing the pace and timing of cortical microstructural development in youth, using submillimeter resolution, B1+ bias-corrected T1w/T2w mapping, and mapping data to a multimodal cortical parcellation in a large, demographically diverse sample.
Spatial pattern of age-related increases in T1w/T2w remain consistent while controlling for cortical thickness
Cortical thickness and myelin content are systematically related, where thinner areas in sensorimotor cortex have higher myelin content, and thicker areas in paralimbic cortex have lower myelin content (but with exceptions such as the thick, heavily myelinated primary motor cortex; Glasser et al., 2014). Previous studies have observed variable rates of cortical thinning across the cortex during early childhood and adolescence (Shaw et al., 2008; Vandekar et al., 2015). Further, age-related changes in cortical thickness have been difficult to distinguish from changes in cortical myelin content under certain circumstances (Natu et al., 2019). In our study, age-related differences in T1w/T2w were not driven by cortical thinning, as we found highly consistent topography of age-related differences in T1w/T2w after controlling for cortical thickness in each parcel. This indicates that cortical thinning and T1w/T2w differences during youth are dissociable developmental processes.
Neurobiological mechanisms: the role of cortical myelin in regulating plasticity
While both myelin-sensitive imaging and postmortem histology have provided evidence for continuous increases in intracortical myelin content into adulthood (R.A. Hill et al., 2018; Hughes et al., 2018; Grydeland et al., 2019), the neurobiological mechanisms promoting cortical myelination throughout the lifespan remain largely unclear. Studies in animal models suggest that cortical myelination provides a mechanism to continuously tune the firing properties of neural circuits to support adaptive behavior within an individual's environment (Makinodan et al., 2012; Mount and Monje, 2017). Axonal myelination enhances neural signaling speed (McDougall et al., 2018) and plays an important role in regulating the timing of sensitive periods in neurocognitive development by providing structural constraints that limit local plasticity (Takesian and Hensch, 2013). Research using animal models has provided strong evidence for the role of myelin-related proteins such as Nogo-A in reducing plasticity by inhibiting neurite outgrowth and dendritic arborization in cortical circuits (Chen et al., 2000; McGee et al., 2005).
A role for cortical myelination in regulating plasticity during early sensitive periods for sensory cortex development has been demonstrated in animal models (Takesian and Hensch, 2013; Toyoizumi et al., 2013). However, the role of cortical myelination in regulating higher-order cognitive development in humans remains poorly understood. Recent work suggests that a large proportion of neocortical myelin ensheathes axons of parvalbumin-positive (PV+) fast-spiking inhibitory neurons (Micheva et al., 2016). Further, the myelination of PV+ interneurons is critical for maintaining excitation-inhibition (E/I) balance in maturing cortical sensory circuits (Benamer et al., 2020). The maturation of PV+ interneurons and E/I balance in neural circuits is a crucial factor regulating the opening and closure of sensitive periods in neurodevelopment (Takesian and Hensch, 2013; Toyoizumi et al., 2013; Takesian et al., 2018). Moreover, we speculate that cortical myelin development could be linked to organizational gradients of PV gene expression and the modulation of inhibitory circuitry throughout the cerebral cortex (Burt et al., 2018).
Limitations
Cortical MR signals are sensitive to many tissue properties including iron, myelin, cell density, and water content, and contrast images such as T1w/T2w generally represent a mix of these properties that can vary nonlinearly across laminar architecture (Carey et al., 2018; Glasser et al., 2022). Further, age-related differences in T1w/T2w could reflect a combination of these cortical tissue changes. For example, in addition to increased intracortical myelin content, age-related increases in T1w/T2w could also reflect an increase in tissue iron concentration (Larsen et al., 2020), the reduction in unmyelinated tissue processes such as dendrites, which could result in higher apparent myelin concentration per unit volume (expected to be associated with cortical thinning), or the dehydration of tissues (as occurs in early development).
Future directions
While this study characterized age-related changes in cortical T1w/T2w during youth aged 8–21 years, the cross-sectional design precluded inferences on intraindividual changes in cortical T1w/T2w. Although the HCP-D study contains a longitudinal component (Somerville et al., 2018), it was not feasible to incorporate the longitudinal measurements in this investigation because of ongoing data collection and processing. A more comprehensive lifespan approach could help identify milestones of cortical myelin development and senescence, expanding on previous work (Grydeland et al., 2019). The current study provides an important reference dataset characterizing adolescent cortical microstructural development that can serve as a benchmark for replication using the longitudinal subset of the HCP-D sample or other developmental longitudinal datasets such as ABCD (Casey et al., 2018). Such efforts would extend our understanding of cortical myelin development to allow us to evaluate the temporal relationships between cortical myelination and improving cognitive abilities. Future work could also investigate how socioeconomic factors, exposure to early life adversity, or other environmental contexts impact the rate of cortical myelination during youth, as well as the multifaceted relationships between developing cortical microstructure, connectivity, and psychological maturation.
Footnotes
This work was supported by National Institute of Mental Health (NIMH) Grants R01MH129493, R24MH108315, R24MH122820, U01MH109589, and U01MH109589-S1; the 14 National Institutes of Health (NIH) institutes and centers that support the NIH Blueprint for Neuroscience Research; the McDonnell Center for Systems Neuroscience at Washington University; the Office of the Provost at Washington University; and a Seed Grant from the Harvard Brain Institute Bipolar Disorder Seed Grant Program. Portions of this research were carried out at the Harvard Center for Brain Science using instrumentation supported by the NIH Shared Instrumentation Grant Program S10OD020039. This work was also supported by NIMH grants R01 MH60974 and the National Institute of Biomedical Imaging and Bioengineering (NIBIB) grant T32 EB021955. We thank Melanie Grad-Freilich for assistance with data collection and data access, Erin Reid for assistance with neuroimaging data quality assurance, and Valerie Sydnor for helpful discussion and sharing data related to the sensorimotor-association axis.
The authors declare no competing financial interests.
- Correspondence should be addressed to Graham L. Baum at gbaum{at}fas.harvard.edu
