Quantitative MRI provides markers of intra-, inter-regional, and age-related differences in young adult cortical microstructure

Measuring the structural composition of the cortex is critical to understanding typical development, yet few investigations in humans have charted markers in vivo that are sensitive to tissue microstructural attributes. Here, we used a well-validated quantitative MR protocol to measure four parameters (R1, MT, R2*, PD*) that differ in their sensitivity to facets of the tissue microstructural environment (R1, MT: myelin, macromolecular content; R2*: myelin, paramagnetic ions, i.e., iron; PD*: free water content). Mapping these parameters across cortical regions in a young adult cohort (18–39 years, N = 93) revealed expected patterns of increased macromolecular content as well as reduced tissue water content in primary and primary adjacent cortical regions. Mapping across cortical depth within regions showed decreased expression of myelin and related processes – but increased tissue water content – when progressing from the grey/white to the grey/pial boundary, in all regions. Charting developmental change in cortical microstructure cross-sectionally, we found that parameters with sensitivity to tissue myelin (R1 & MT) showed linear increases with age across frontal and parietal cortex (change 0.5–1.0% per year). Overlap of robust age effects for both parameters emerged in left inferior frontal, right parietal and bilateral pre-central regions. Our findings afford an improved understanding of ontogeny in early adulthood and offer normative quantitative MR data for inter- and intra-cortical composition, which may be used as benchmarks in further studies.

Measuring the structural composition of the cortex is critical to understanding typical development, yet few investigations in humans have charted markers in vivo that are sensitive to tissue microstructural attributes. Here, we used a well-validated quantitative MR protocol to measure four parameters (R 1 , MT, R 2 *, PD*) that differ in their sensitivity to facets of the tissue microstructural environment (R 1 , MT: myelin, macromolecular content; R 2 *: myelin, paramagnetic ions, i.e., iron; PD*: free water content). Mapping these parameters across cortical regions in a young adult cohort (18-39 years, N ¼ 93) revealed expected patterns of increased macromolecular content as well as reduced tissue water content in primary and primary adjacent cortical regions. Mapping across cortical depth within regions showed decreased expression of myelin and related processesbut increased tissue water contentwhen progressing from the grey/white to the grey/pial boundary, in all regions. Charting developmental change in cortical microstructure cross-sectionally, we found that parameters with sensitivity to tissue myelin (R 1 & MT) showed linear increases with age across frontal and parietal cortex (change 0.5-1.0% per year). Overlap of robust age effects for both parameters emerged in left inferior frontal, right parietal and bilateral precentral regions. Our findings afford an improved understanding of ontogeny in early adulthood and offer normative quantitative MR data for inter-and intra-cortical composition, which may be used as benchmarks in further studies.
A core challenge for human neuroscience is the design of robust anatomical imaging methods that are sensitive to inter-regional differences in tissue properties, and to profiles of intra-cortical tissue change from the grey-white border to the pial surface in any one region. The parcellation of human cortex based on cyto-and myeloarchitectonic boundaries has been a major pursuit since the work of Brodmann and Flechsig in the early 20th century (Sereno et al., 2013;Glasser et al., 2016;Turner, 2015;Nieuwenhuys, 2013;Nieuwenhuys et al., 2014;Zilles et al., 2015). However, it is only recently that such questions have been addressed in-vivo in humans. This is made possible by the use of magnetic resonance imaging (MRI), which can provide data for morphometry (Ashburner and Friston, 2000;Dale et al., 1999;Fischl et al., 1999aFischl et al., , 1999b or microstructure (Weiskopf et al., 2015).
The MR signal is sensitive to many important tissue properties, such as iron content, myelin, cell density and water content; however, the contrast-weighted images (T 1w , T 2w ) typically used in MRI reflect a complex mix of these properties that can vary non-linearly across the imaged volume. By comparison, Quantitative MRI (Bock et al., 2013;Barazany and Assaf, 2012;Dinse et al., 2013Dinse et al., , 2015Stüber et al., 2014;Marques et al., 2010; for review, see Turner, 2015Turner, , 2016Bazin et al., 2014;Cohen-Adad, 2014;Sereno et al., 2013;Dick et al., 2012) can be used to map specific MRI properties of tissue in order to provide indices of microstructure, myelination and related cellular processes (Helms et al., 2008aWeiskopf et al., 2013;Lutti et al., 2014) in a timeefficient manner with high spatial specificity. It thus provides the opportunity to acquire a multi-modal, whole-brain view of developmental changes in underlying tissue properties.
Charting the normal development and aging of human cortical tissue is a fundamental goal of neurobiology, and is also critical for accurately characterizing atypical development, individual differences, and shortand long-term plasticity. Development is reported to follow a posteriorto-anterior gradient with primary areas maturing earliest in life and association areas, which mediate higher-order functions, developing later (Gogtay et al., 2004; for review of processes, see Marsh et al., 2008). At earlier points in development through adolescence, there is evidence to suggest that deviation from typical trajectories may increase vulnerability to psychiatric disorders (Thompson et al., 2001;Greenstein et al., 2006;Sowell et al., 2003;Shaw et al., 2007) whereas in later life, such deviation may be indicative of neurodegenerative decline, for which age is often the greatest predictor (Barkhof et al., 2009;Bartzokis, 2004Bartzokis, , 2011Frisoni et al., 2010). A cornerstone in the development of mature cortex is the emergence of myelinated fibres within the cortical sheet (Flechsig, 1920;Yakovlev and Lecours, 1967;Deoni et al., 2015). Though the exact trajectories are unclear, the rate at which change occursand the age at which development stabilizesare thought to be region-specific (e.g. Yeatman et al., 2014;Whitaker et al., 2016) and to interact with functional organization (e.g. Yeatman et al., 2012;Gomez et al., 2017).
To date, few quantitative imaging studies have explored developmental changes in tissue composition across cortex from late adolescence to the mid-thirties. This is a crucial age range to characterize, not least because it is the 'sample of choice' for the vast majority of structural and functional MRI studies. Here, we used the multi-parameter mapping (MPM) protocol Callaghan et al., 2014aCallaghan et al., , 2014bHelms et al., 2008aHelms et al., , 2008b to explore potential parameter-specific (R 1 , MT, PD*, R 2 *) variation in tissue over the depth of the cortical sheet, and across a range of cortical regions. Further, we charted age-related differences in cortical microarchitecture across early adulthood. We mapped a set of normative, cortical-depth-specific regional MPM values for young adults that can be used as reference values for future studies. Moreover, we found considerable, region-specific age-related changes in parameters related to the degree of tissue myelin content and myelinrelated processes.

Participants
Participants were 93 right-handed healthy adults (mean age ± SD: 23.6 ± 4.3; range: 18-39; 57 female, 36 male). The study received approval from the local ethics committee. All scanning took place at the Wellcome Trust Centre for Neuroimaging (WTCN), London.
Participants were sampled over approximately 24 months. Thirtyfour participants were recruited as part of a study of musicianship and consisted of expert violinists (n ¼ 18; mean age ± SD: 22.8 ± 2.8; 13 female, 5 male) and closely matched non-musicians (n ¼ 16; mean age ± SD: 23.3 ± 3.1; 12 female, 4 male). All had completed or were enrolled in a university degree, and were recruited from the University of London, music conservatories in London, and local participant pools. We analyzed data for effects of violin expertise and will report these findings in a subsequent report. In brief, effects of violin expertise in cortex were modest and emerged only in ROI analyses of primary auditory cortex, where we found limited evidence of significant age-related effects in the present study.
The remaining participants (n ¼ 59; mean age ± SD: 23.9 ± 4.9; 32 female, 27 male) were sampled from the general population through local participant pools. These subjects took part in three experiments: one exploring the potential association between auditory perceptual abilities, musicianship, tonotopic organization and structural properties of the auditory cortex (data not reported here), one investigating the relationship of trait empathy and brain microstructure (Allen, Frank, et al., 2017), and a third investigating metacognition and MPM assays (Allen et al., 2017).
There was no significant difference in age between genders across the full cohort (z ¼ 0.85, p > 0.4), nor any significant effects of gender on MPMs in any models (all p > 0.3).

Data acquisition
The multi-parameter mapping protocol data Lutti et al., 2010Lutti et al., , 2012 were acquired at the WTCN using a 3T wholebody Tim Trio system (Siemens Healthcare) with radiofrequency body coil for transmission and a 32-channel head coil for signal reception. The MPM protocol consisted of three differently weighted 3D multiecho FLASH acquisitions acquired with 800 μm isotropic resolution. Volumes were acquired with magnetization transfer (MT w ), T 1 -(T 1w ), and proton density (PD w ) weighting. The MT weighting was achieved through application of a Gaussian RF pulse (4 ms duration, 220 nominal flip angle) applied 2 kHz off-resonance prior to non-selective excitation.
The sequence settings of the MPM protocol were modified following collection of data for the musicianship sample (cohort 1, n ¼ 34), reflecting the on-going development of the MPM sequences at the WTCN. Cohort-specific details follow.

Cohort 1
A field of view of 256 Â 224 Â 166 mm 3 was used with a matrix size of 320 Â 280 x 224. The PD w and T 1w volumes were acquired with a TR of 25.25 ms and a flip angle of 5 and 29 respectively. The MT w volume was acquired with a TR of 29.25 ms and excitation flip angle of 9 . The excitation employed a hard pulse; RF spoiling was used with an increment of 50 , and gradient spoiling producing 6*pi dephasing across a voxel. The acquisition was accelerated by using GRAPPA (with a parallel imaging factor of 2 with 18 integrated reference lines) in the phaseencoded direction (AP) and by a partial Fourier acquisition in the partition direction (RL, with factor 6/8). To improve image quality (maximize SNR and minimize geometric distortion at the same time), eight gradient echoes with alternating readout polarity were acquired with high readout bandwidth (460 Hz/pixel) with echo times ranging from 2.39 ms to 18.91 ms in steps of 2.36 ms.
A fixed modification of MPM image slab orientations (30 ) was applied for some subjects to counter image artifact due to eye movement and blinking; this change in the acquisition did not yield any significant differences (testing binary main effect of slab rotation present/absent; p < 0.01, whole-brain uncorrected) for any map (R 1 , PD, MT, R 2 *) between subjects with and without slab rotation (see footnote 3). Image slabs for field maps were all non-rotated along the axial orientation.
The B 1 mapping acquisition consisted of 15 measurements with nominal flip angle ranging from 135 to 65 in 5 decrements. The total scanning time of the MPM protocol was approximately 37 min.

Cohort 2
For the second cohort, the MPM protocol was modified to improve accuracy by accounting for non-linearities in the transmit chain . To achieve this, different flip angles for the PD-, MT-(both 6 ) and the T 1 -(21 ) weighted acquisitions were achieved by scaling the duration of the pulse while maintaining a constant B 1þ amplitude (i.e. a consistent operating point for the RF amplifier) that additionally matched that used for the B 1 þ mapping sequence. Gradient echoes were again acquired with alternating readout gradient polarity using a readout bandwidth of 488 Hz/pixel. Eight equidistant echo times ranging from 2.34 to 18.44 ms in steps of 2.3 ms were acquired for the PD w and T 1w acquisitions. Only the first six echoes were acquired for the MT-weighted acquisition in order to maintain a 25 ms TR for all of the FLASH volumes. To further accelerate the data acquisition, the partial Fourier acquisition scheme in the partition direction was replaced by parallel imaging with an acceleration factor of 2 again using the GRAPPA algorithm, now with 40 integrated reference lines in each phase-encoded direction. A 30 slab rotation was used for all acquisitions in this cohort.
The B 1 mapping acquisition consisted of 11 measurements with nominal flip angle ranging from 115 to 65 in 5 decrements. The total scanning time of the MPM protocol was approximately 26 min.

Procedure
Participants provided written informed consent and were screened for contraindications for MRI. B 1 þ and B 0 field maps were collected at the beginning of each session, followed by the MT, PD w , and T 1w scans. Participants' eye and head movements were monitored using an eye tracker (Eyelink 1000 Core System) during scanning runs. Rest breaks of several minutes were provided between scans as required.
Data pre-processing Images were pre-processed using the Voxel Based Quantification (VBQ) toolbox in SPM 8. In brief, regression of the log signal from the echoes of all weighted volumes were used to calculate a map of R 2 * using the ordinary least squares ESTATICS approach . The set of echoes for each acquired weighting were then averaged to increase the signal-to-noise ratio (Helms and Dechent, 2009). This was done using only the first six echoes for Cohort 2. The 3 resulting volumes were used to calculate MT, R 1 , and PD* maps as described in Helms et al. (2008aHelms et al. ( , 2008b and Weiskopf et al. (2013). Quantitative R 1 values at each voxel were estimated based on the rational approximation of the Ernst equation described by Helms et al. (2008a). To maximize the accuracy of the R 1 map, these maps were corrected for transmit field inhomogeneities by constructing a map from the calibration data according to the procedure detailed in Lutti et al. (2012). The R 1 maps were also corrected for imperfect spoiling characteristics using the approach described by Preibisch and Deichmann (2009). The MT map was constructed using the procedure described in Helms et al. (2008b). This is a semi-quantitative metric depicting the percentage loss of magnetization resulting from the MT pre-pulse used and differs from the commonly used MT ratio (percentage reduction in steady state signal) by explicitly accounting for spatially varying T 1 relaxation times and flip angles . Finally, PD* maps were estimated from the signal amplitude maps by adjusting for receive sensitivity differences using a post-processing method similar to UNICORT (Weiskopf et al., 2011). To make the PD*maps comparable across participants, they were scaled to ensure that the mean white matter PD* for each subject agreed with the published level of 69% (Tofts, 2003). This quantity is referred to as effective PD (PD*) because it was calculated based on the average FLASH volumes and there was no correction for R 2 * signal decay.
Following reconstruction of multi-parameter images, all images were manually inspected for any evidence of alignment difficulties, head movement or other image artifacts (e.g., aliasing) by a rater who was blind to subject identity.

Cortical surface reconstruction
Participants' cortical surfaces were reconstructed using FreeSurfer (v. 5.3; Dale et al., 1999). Use of multi-parameter maps as input to Free-Surfer can lead to localized tissue segmentation failures due to boundaries between the pial surface, dura matter and CSF showing different contrast compared to that assumed within FreeSurfer algorithms . Therefore, an in-house FreeSurfer surface reconstruction procedure was developed to overcome these issues. Full details of the processing pipeline are provided in supplemental methods.

MPM data extraction
First, all subjects were rotated to the same (canonical) orientation, using the AFNI 3dwarp routine (-deoblique flag). MPM data were then 3 For the first MPM cohort (n ¼ 34), T 1w , PD w and MT images were acquired with different slab orientations over some subjects relative to others. Initial inspection of data acquired with the slab aligned to each cardinal axis showed susceptibility artifact that affected cortex in a subset of participants. Although eye movements were monitored during scanning runs, slight movement (e.g., due to blinking) led to artifact within orbitofrontal and medial temporal lobes in some datasets. To counter this issue, the acquisition protocol was modified, by rotating each image slab at 30 about the x-axis (such that the eyes lay outside the slab). Participants with data acquired without slab rotation were inspected blind to subject identity for evidence of susceptibility artifact; those participants that showed evidence of artifact within cortical areas were re-scanned using the rotated acquisition protocol. In total, 6 of 34 participants showed susceptibility artifact with the unrotated acquisition and were re-scanned with the rotated protocol; 15 of 34 participants showed no evidence of susceptibility artifact with the original unrotated acquisition and were not re-scanned; 13 of 34 participants were scanned using the rotated protocol as default. Whole-brain analyses of each MPM using slab rotation as a regressor of interest showed no significant differences across any vertices over either hemisphere for participants with rotated versus unrotated acquisition (p < 0.01, uncorrected). Moreover, there was no significant difference in age between the subjects with non-rotated versus rotated acquisition (z < 0.01, p > 0.99). mapped onto each subject's surface (using the FreeSurfer mri_vol2surf routine). For each reconstructed hemisphere, quantitative data were sampled along the normal to each surface vertex, for cortical depth fractions from 0.1 (i.e., above white matter surface boundary) to 0.9 (i.e., beneath pial surface boundary) in increments of 0.1 (see Dick et al., 2012).

Analyses
For each relaxation parameter (R 1 , MT, R 2 *, PD*), we first created cross-subject hemisphere-wise average maps for each cortical depth sampling fraction (0.1-0.9) using cortical-surface-based methods with curvature-based alignment (Fischl et al., 1999a,b;Dick et al., 2012;Sereno et al., 2013). In each of a series of cortical regions-of-interest, selected to sample from a range of primary, secondary, and tertiary areas across the brain, we extracted hemisphere-wise mean estimates of each relaxation parameter across depth fractions. Cortical regions were defined from a standard FreeSurfer atlas (aparc.2009), and sampled onto each subject's cortical surface during reconstruction. Regions were: superior pre-central sulcus; subcentral gyrus/sulcus; inferior pre-central sulcus; inferior frontal sulcus; middle frontal sulcus; superior parietal gyrus; angular gyrus; Heschl's gyrus; inferior occipital gyrus/sulcus; superior temporal gyrus/planum temporale; probabilistic area MT; probabilistic V1; posterior collateral sulcus; superior occipital gyrus; parieto-occipital sulcus; subparietal sulcus; middle cingulate gyrus/sulcus. ROI mean estimates for each relaxation parameter were averaged across hemispheres at each cortical depth fraction sampled (0.1-0.9). Previous histological and MR literature has demonstrated non-linear change in tissue properties related to myelin processes, progressing from the grey-white to grey-CSF boundary (Annese et al., 2004;Walters et al., 2003;Whitaker et al., 2016;Stüber et al., 2014;Waehnert et al., 2016;Dinse et al., 2015). Here, multi-level models were fitted to ROI mean estimates for each MPM to model linearity versus non-linearity of parameter change across depth fractions, whilst accounting for variance across subjects and ROIs. For each parameter, we specified fixed effects of MPM cohort, ROI, gender, as well as interaction terms for depth that specified quadratic or cubic fits. Each model also included random effects of subject and ROI (i.e., three-level hierarchy: depth fractions nested within ROIs, and ROIs nested within subjects). We initially tested models with and without linear random (i.e., personspecific) slopes across depth fractions (i.e., nested within ROIs), and found that fits were always significantly improved with the addition of linear random slopes (all p < 0.0001, likelihood ratio test); thus linear random slopes were specified in all models. Models were first specified for each MPM with linear fixed effect terms for depth, and latterly with quadratic and cubic interaction fixed effect terms for depth. Likelihood ratio tests were used to compare model log-likelihoods, testing whether addition of each non-linear term yielded any robust change in model fit, relative to the simpler alternative model (i.e., linear vs. quadratic, quadratic vs. cubic). All multi-level models were fitted in STATA 14 (STATA Corp).
For age-based analyses, we used participants' age in whole years as a continuous linear regressor at each vertex per hemisphere, with initial analyses carried out in Qdec (FreeSurfer v. 5.3). For age analyses, each subject's data were smoothed with a surface (2D) kernel of 10 steps (approximating a Gaussian 2D kernel of 3 mm FWHM; . Vertex-wise age analyses are reported sampling at 0.5 cortical depth. To avoid inflation of type-1 error that would have resulted from running vertex-wise analyses for all cortical depth fractions, this depth fraction was selected as the most representative of mid-cortical depth profiles (Dick et al., 2012;Sereno et al., 2013;Lutti et al., 2014;cf. Waehnert et al., 2014). However, we also explored potential interactions of age effects across depth fractions and ROIs (further to Whitaker et al., 2016) using multi-level models. Specifically, the following interaction terms were added to the models described above as fixed effects: age x depth; age x ROI; age x depth x ROI.
Previous studies have found that R 1 values are associated with variation in the local curvature and thickness of the cortex; thus, R 1 measurements tend to be increased in thicker, more highly convex regions (e.g., gyral crowns) (Sereno et al., 2013;Dick et al., 2012;Waehnert et al., 2014). In addition, changes in the MPM acquisition protocol between the cohorts we scanned here (see 2.2.1 & 2.2.2) were associated with differences in MPM map values. The second cohort, scanned with the protocol that better addressed non-linearities in the transmit chain, showed consistently greater (up to~15%) R 1 values with a gently spatially varying pattern that could be reproduced by comparing results in a single individual subject scanned across both protocols. Therefore, to control for these effects, we regressed out local curvature, cortical thickness and MPM cohort (see also Grydeland et al., 2013); we then performed Pearson correlations at each vertex between age and MPM values that were residualized by curvature, thickness and MPM-cohort. In general, curvature-, thickness-, and cohort-residualized age regressions were similar to age regressions using raw parameter values.
To produce unbiased estimates of age effects on each parameter we used a 'leave-one-out' jackknife procedure. Here, we first performed vertex-wise age-MPM Pearson correlations for the full cohort, and then repeated the procedure iteratively omitting one subject in each instance. Pearson r-values for the full cohort and each leave-one-out partial estimate were Fisher z-transformed, and a mean of partial estimates calculated. Jackknife estimates at each vertex were calculated as: (N)(T) -(N-1)(Tm), where N ¼ 93; T ¼ full cohort z-transformed rvalue; Tm ¼ mean of partial estimates (z-transformed before averaging). Finally, vertex-wise Jackknife estimates were re-transformed to r-values and corresponding p values were calculated. The jackknifing procedure was performed in regression models with age as a vertexwise predictor of raw MPM values, and also in regression models with age as a vertex-wise predictor of cohort-, thickness-and curvatureresidualized maps.
Jackknifed statistical maps were thresholded using peak-level False-Discovery Rate (FDR) correction (Benjamini and Hochberg, 1995); FDR-corrected q < 0.05, per hemisphere. For illustrative purposes, we also identified regions where significant jackknifed effects of age overlapped for both the R 1 and MT multi-parameter maps. Note that these maps show greatest sensitivity to cortical myelin, and thus were of central interest here (see Lutti et al., 2014;Sereno et al., 2013). Per hemisphere, we determined the vertices that survived FDR-correction (q < 0.05) for jackknifed age analyses of the cohort-, thickness-and curvature-residualized MPMs (i.e., R 1 & MT). Using Matlab, we then created a binary mask per hemisphere corresponding with vertices where the jackknifed model results for both MPMs showed FDR-significant effects of age. Clusters of vertices reflecting the unsmoothed overlap of the age effects for the two MPMs (R 1 & MT) were extracted and defined as ROIs on a standard cortical surface; ROIs were sampled onto each subject's cortical surface. Across each of these ROIs, we plotted the linear relationship between age and subject-wise ROI means for R 1 and MT.
Cortical-surface-averaged group data for each parameter as well as all individual subject data and analysis scripts are available at: https://doi. org/10.18743/DATA.00011.

Results
Here, we explored intra-and inter-regional differences in quantitative markers of tissue microstructure, using a multi-parameter mapping protocol that affords quantitative MR proxies for myelin-related tissue processes. Further, we used a cross-sectional design to explore the effects of age across early adulthood on cortical myelination. To develop a clearer understanding of inter-and intra-regional myelin-related tissue properties, we produced average maps for each of the MPMs across the cortical sheet, sampling halfway through cortex; moreover, we charted differences in MPM parameters over the depth of the cortical sheet, within and between a range of cortical ROIs.

Group average MPM results
Average MPMs (Fig. 1) revealed the expected inter-regional differences in cortical myelin and myelin-related processes, in line with previous literature. Parameters that show greatest sensitivity to myelin and related processes (R 1 , MT, and R 2 *), had highest values within motor and sensory regions, including somatomotor, auditory, and visual cortex. Further to previous studies (Glasser and Van Essen, 2011;Glasser et al., 2016;Sereno et al., 2013;Waehnert et al., 2016;Bock et al., 2013;Mangeat et al., 2015), R 1 , MT and R 2 * revealed strips of dense myelination over pre-central and post-central regions, presumably in areas 4 and 3b/1, respectively, with an intervening lowermyelin septum (likely area 3a; Geyer, 2013;Dinse et al., 2015;Stüber et al., 2014; see also Glasser and Van Essen, 2011). R 1 , MT and R 2 * also revealed the heavily myelinated auditory core at the most medial aspect of Heschl's gyrus (Dick et al., 2012;Sigalovsky et al., 2006;de Martino et al., 2015); planum temporale and parts of lateral superior temporal gyrus (STG) additionally showed elevated R 1 and MT values (see Glasser and Van Essen, 2011;Glasser et al., 2016;Sigalovsky et al., 2006). In visual areas, R 1 and R 2 * exposed the heavily myelinated V1 extending across the calcarine sulcus (Sereno et al., 2013;Fracasso et al., 2016;Cohen-Adad et al., 2012); however, parameter MT showed high values that were restricted to gyral banks flanking the calcarine sulcus (Fig. 1, parameter MT medial surface panels). A possible source of this difference is the local contribution of myelin differences to macromolecular effects at gyri, versus the more anatomically diffuse effects of iron associated with oligodendrocyte cell bodiescontributing in part to the R 1 signal, via R 2 *found across sulci (e.g., Stüber et al., 2014; see 3.3, and discussion, 4.1).
Alternatively, local patterns of cortical folding and curvature may have influenced the detection of macromolecular content by MT in highly concave cortical regions; cortical thickness in sulcal depths (and indeed much of V1) is roughly twice the voxel dimensions, and therefore the contribution of the thin myelinated layers to the overall contrast will be attenuated.
R 1 , MT and R 2 * additionally revealed higher visual areas including V6 (medial surface, dorsal to V1), V3/V3a (lateral surface, dorsal to V2), and area MT (proximal to postero-lateral bounds of inferior temporal sulcus). R 1 and MT also revealed several heavily myelinated cortical regions posterior to post-central gyrus, likely including multi-modal VIP and LIP (Sereno et al., 2013;Glasser and Van Essen, 2011;Waehnert et al., 2016;Huang et al., 2012). Foci of high R 1 and MT potentially co-located with the frontal eye fields also emerged, lying proximal to the dorsal-most aspect of the middle frontal gyrus and the boundary with pre-central gyrus (Glasser et al., 2016;Glasser and Van Essen, 2011;Sereno et al., 2013). PD* revealed a more distinct patterning of regions than the other MPMs, reflecting its high affinity for unbound protons (i.e., tissue water; Baudrexel et al., 2016). Regions typically low in myelination (see Nieuwenhuys, 2013;Nieuwenhuys et al., 2014;Geyer, 2013;Zilles et al., 2015) tended to show highest PD* values, including: the circular sulcus and adjoining insular cortex; cingulate gyrus and sulcus; medial prefrontal cortex; collateral sulcus; parieto-occipital sulcus and regions anterior to V1; superior temporal sulcus; middle temporal gyrus; inferior frontal sulcus; and presumptive Area 3a. PD* also tended to show highest values in sulci and locally concave cortex, a rough mirror image of the distribution for the other three parameters, and likely due to the lower overall myelination that has been noted in sulci.

Intra-regional MPM results
Across depth fractions, we observed the expected pattern of decrease in R 1 and MT values over all cortical ROIs ( Fig. 2a and b). Reduction in R 1 and MT values broadly followed a decaying pattern, with sharpest decreases at depth fractions close to the white matter and pial surfaces (0.1 and 0.9, respectively; see Fig. 2a and b). In line with previous work (Annese et al., 2004;Dinse et al., 2015), this likely reflects the reduction in myelinated fibre density with progression away from the white matter surface toward the pial surface. In modeling the reduction in R 1 and MT values over the cortical sheet in each ROI, we found that cubic trends provided the best fit in every region, versus linear and quadratic trends (p < 0.0001, likelihood ratio test vs. quadratic model).
R 2 * values also decreased across all ROIs progressing from the white matter to the pial surface. Modeling the reduction in R 2 * values, quadratic trends provided the best fit to data across depth fractions in all ROIs (versus linear and cubic trendsboth p < 0.0001). In many ROIs, R 2 * values showed an upward trend after the initial decrease, particularly Panels present mean ± 95% CI (CI across individuals), for each parameter across depth fractions within each ROI (ROI means averaged across hemispheres), extending from proximal to grey-white boundary (gm/wm; 0.1) to proximal to grey-pial boundary (gm/pial; 0.9) (x axes, left to right). Line color denotes the ROI (see inset: corresponding ROI color displayed on standard inflated cortical surface; color-coded ROI names listed at centre). (b) Depth profile plots for medial and ventral surface ROIs. All other specifications as per (a). at cortical depths close to the pial surface (i.e., 0.8 and 0.9 cortical depth). This pattern was noted in regions including probabilistic MT/V5, inferior occipital gyrus/sulcus, inferior frontal and inferior pre-central sulci, and angular gyrus (see Fig. 2a), as well as probabilistic V1, superior occipital gyrus, and subparietal sulcus (see Fig. 2b). A likely cause of this is the presence of blood vessels close to the pial surface, where high levels of iron within hemoglobin would elevate observed R 2 * values (see footnote 4).
PD* values showed increases toward the pial surface from the white matter surface, across all ROIs. Similar to R 1 and MT values, change in PD* values was greatest at depth fractions proximal to the white matter and pial surfaces (0.1 and 0.9, respectively), with PD* values showing high growth there (cf. the declines noted in R 1 and MT values). In all ROIs, cubic trends provided the best fit to PD* increases over depth fractions (fits specified as above for R 1 and MT).

Inter-regional MPM results
Across cortical ROIs, we found clear differences in the mean values of R 1 , MT, R 2 * and PD* parameters, marked by substantial changes in those parameters when sampling across primary and non-primary cortical ROIs.
PD* tended to be reduced in heavily-myelinated areas (e.g., Heschl's gyrus; Fig. 2a, mauve trace). However, we note that inter-regional PD* curves did not reflect a strict 'mirror-image' of areas with elevated R 1 and MT values. For instance, while superior pre-central sulcus (Fig. 2a, light green traces) showed elevated R 1 and MT curve values, PD* curve values in this region were also increased compared to adjacent cortical areas (see Fig. 2a, PD* panel).
Age effects on R 1 /MT/R 2 * Exploring development of myelin and related tissue processes crosssectionally, we correlated age in years with vertex-wise R 1 , MT and R 2 * values that had been residualized by cortical thickness, curvature and MPM cohort (see 2.4.2). To limit bias in model fits, we calculated 'leaveone-out' estimates of models per MPM using a jackknifing procedure (see 2.4.2).
We found evidence of significant (FDR q < 0.05 per hemisphere) correlations between age and R 1 , and age and MT, across the lateral cortical surface. However, we did not find any regions where age and R 2 * correlations survived with FDR-correction (q > 0.05) (but see below).
Positive age-R 1 correlations (i.e., increasing R 1 with age; Fig. 3a) were widespread, extending across much of pre-frontal, frontal and parietal cortex bilaterally. These positive age-R 1 correlations were observed bilaterally at pars opercularis, middle frontal gyrus (MFG), pre-central gyrus and central sulcus, superior frontal gyrus, and superior and inferior parietal lobules (including angular and supramarginal gyri). Additional positive correlations were found in right occipito-temporal regions, including a peak proximal to visual area MT/V5. A further peak emerged at right superior temporal sulcus.
Significant positive age-MT correlations were less extensive than those observed for age-R 1 (see Fig. 3b). The largest age-MT peaks manifested at left pre-central gyrus and central sulcus, along with a series of peaks across right pre-central gyrus. Other smaller peaks emerged at right supramarginal gyrus, left pars opercularis, right MFG, and right superior parietal lobule.
Exploring age effects within ROIs, we fitted multi-level models testing fixed main effects of MPM cohort, gender, age, ROI and depth fraction (depth fitted with linear and non-linear terms), along with age x ROI and age x ROI x depth fraction interactions. MT, R 1 and R 2 * each showed robust main effects of the age and depth fraction terms, with significant differences found amongst ROIs (see Supplemental Tables 1-3). Notably, although age effects on R 2 * were not robust in the vertex-wise analysis (see above), age was a significant predictor of R 2 * in the multi-level analysis (slope: 0.053 [~0.3%/year]). Importantly, age Â depth interactions for each parameter were also significant (all p < 0.01); in line with previous studies (Whitaker et al., 2016), this suggests that effects of age were attenuated moving closer to the pial surface. Age x ROI Â depth interactions also reached significance for each parameter; however, the coefficients associated with these three-way interactions were small, suggesting the change in age effects across depth fractions and between ROIs were subtle. Supplemental Tables 1-3 summarize these effects for MT, R 1 and R 2 *; effects are plotted across age quintiles in Supplemental Fig. 1.
Finally, Supplemental Fig. 3 shows a map of age-related changes in cortical thickness (all decreases) with cohort and slab rotation as nuisance factors; along the inferior, middle, and superior frontal gyri as well as the supramarginal gyrus, these changes were concordant in part with those reported by Whitaker et al. (2016) and Hogstrom et al. (2012) for younger and older cohorts with large age ranges, but unlike these studies did not show any age-related thickness decreases posteriorly. 4 To our knowledge, ours are the first in vivo data to demonstrate elevation of R 2 * values close to the pial surface. A number of studies have explored mapping of T 2 */R 2 * data histologically (e.g., Stuber et al., 2014) or in vivo (e.g., Cohen-Adad et al., 2012;Marques et al., 2017). Many studies that have charted R 2 * measurements across cortex have done so sampling at a single fixed cortical depth (Cohen-Adad et al., 2012), or at only a small number of points through cortex . Marques et al., 2017 did observe a slight plateauing of R2* at 20% of cortical depth (toward the pial surface) in some ROIs at 7T. Regarding comparisons of post-mortem histology/MR to in vivo imaging, a challenge is the removal of vasculature from the brain in the ex vivo state. This means that any direct comparison of measurements from ex vivo MR or histology to those made with MR in vivo, must be viewed in light of the removal of a potential source of signal variation (i.e., major vessels) in ex vivo analyses. It may be possible in future analyses to construct maps of vasculature (for instance, using PD* data), and in turn to mask these out from ROIs, as a means of appraising the impact of vessels on R 2 * values near the pial surface.
D. Carey et al. NeuroImage 182 (2018) 429-440 Overlap of age effects -R 1 and MT As a way to explore the extent to which age effects were common to both R 1 and MT, we compared age effects for overlapping significant vertices in both analyses. We isolated vertices over each hemisphere where both R 1 and MT (corrected for covariates) had shown FDRsignificant jackknifed age correlations (i.e., vertex-wise q < 0.05), and then defined these regions as ROIs, which we sampled onto each subject's surface (see 2.4.2). Fig. 4 presents age (range: 18-39 years) regressed against the ROI mean R 1 or MT value per subject, over each hemisphere. (Note that the definition of overlap here is quite conservative in that no smoothing was applied to the FDR-corrected cross-subject correlation maps; therefore there will be greater apparent overlap in the surfacesmoothed maps (10 steps) shown in Fig. 3).
Across the left hemisphere (Fig. 4a), we found three regions of overlap for age-R1 and age-MT effects; these encompassed pars opercularis (green), lateral pre-central gyrus/central sulcus (cyan), and ventral central sulcus (blue). Linear regression fits for R 1 showed that across these ROIs, R 1 values increased at a rate of 0.003 s À1 per year. Similarly, over these ROIs, MT values increased at rates ranging from 0.006 to 0.009 pu per year. Across the right hemisphere (Fig. 4b), we identified four regions where age-R1 and age-MT effects tightly overlapped. These included: dorsal (yellow), lateral (purple), and ventral (orange) pre-central gyrus, and supramarginal gyrus (red). Similar to the left hemisphere, R 1 values in the right hemisphere ROIs increased at rates of 0.003-0.004 s À1 per year; MT values increased at rates ranging from 0.006 to 0.008 pu per year.

Discussion
A fundamental challenge in human neuroscience is the development of efficient and robust anatomical imaging techniques that can enable specific tissue properties within the cortex to be quantified in vivo. Such approaches are critical to charting healthy human brain structure across development, together with informing understanding of tissue deterioration in aging and disease (Yeatman et al., 2014;Deoni et al., 2015, Fig. 3. Results of age-MPM correlations. Maps present vertices where positive age-parameter correlations (jackknifed) were significant, after adjusting for effects of local cortical curvature and thickness, and MPM cohort (see Materials and Methods). (a) Significant age-R 1 correlations emerged across much of the lateral surface, including frontal, parietal and temporal regions. (b) Significant age-MT correlations were less extensive, emerging at pre-central, dorso-lateral pre-frontal, and parietal regions. Heat scale overlays indicate range of jackknifed Pearson r coefficients, thresholded by FDR-corrected significance of age-parameter correlations (all effects hemisphere-wise FDR-corrected, q < 0.05). Pearson r maps are displayed with surface smoothing of 10 steps (approximating a 2D Gaussian kernel of 3 mm FWHM). Fig. 4. Overlap of R 1 -and MT-age effects over (a) left and (b) right hemispheres. Insets present ROIs, defined as regions where overlap of age-R 1 and age-MT correlations manifested (at vertex-wise FDR-corrected significance when adjusting for effects of local cortical curvature and thickness, and MPM cohort). Overlap was determined based on comparing vertex-wise FDR-robust effects per parameter, and hence captures the most robust peaks from Fig. 3 (Fig. 3 presents results with surface smoothing of 10 steps). Surrounding panels display scatter plots of age and subject-wise ROI mean parameter values (ROI correspondence denoted by color coding), with age-parameter linear fits (dashed traces; intercept and age slope from linear fits reported). Crosses: R 1 ; circles: MT. 2011; Callaghan et al., 2014a;Laule et al., 2004). Here, we used a previously well-validated quantitative multi-parameter mapping (MPM) protocol (Helms et al., 2008a(Helms et al., , 2008bWeiskopf et al., 2011Weiskopf et al., , 2013Lutti et al., 2010;Callaghan et al., 2014b) to characterize tissue profiles in-vivo across a range of cortical regions. Moreover, we explored effects of aging across late adolescence and early adulthood on myelination, using MPM indices that offer tissue-specific proxies for cortical myelin processes. We found that across a series of cortical ROIs, MPMs sensitive to myelin (R 1 , MT) showed enhanced values in line with expected differences in myelination between primary and primary adjacent regions, as compared to many non-primary cortical regions. Further, we found that depth profiles of MPMs reflective of myelin and myelin-related processes (R 1 , MT, R 2 *) showed a monotonic decline in parameter values across the cortex, when progressing from the white matter surface towards the pial surface. In contrast, the MPM reflective of tissue water (PD*) increased as a function of distance from the white matter surface. Cross-sectional effects of age on myelination were most robust for R 1 , with linear fit against age from vertex-wise analyses showing R 1 increases of~0.5% per year across much of pre-frontal, frontal and parietal cortex. Linear fits of MT against age from vertex-wise analyses estimated an increase of~1.0% per year, but with correlations that were somewhat less extensive spatially than those found for R 1 .
Inter-regional MPM normative data and intra-regional depth profiles A central aim of our study was to use MPMs to explore myelin and myelin-related variation in cortex using MPMs in order to provide normative mapping data from a healthy sample. In line with previous histological (Annese et al., 2004) and combined histological/MR investigations (Walters et al., 2003;Stüber et al., 2014;Fracasso et al., 2016), we found that MPMs allowed us to distinguish between cortical regions (e.g., primary vs. association cortex), and also to chart intra-regional tissue properties, based on profiles of tissue-sensitive parameters through the cortical sheet.
Previous investigations have explored inter-regional and cortical depth profiles for R 1 (¼1/T 1 ) alone, which has typically shown high sensitivity to heavily-myelinated cortical tissue. In particular, Sigalovsky et al. (2006) charted cortical R 1 values across anatomical subdivisions of Heschl's gyrus, while Dick et al. (2012) characterized the cortical depth profile of R 1 values within sub-parcellations of Heschl's gyrus (Te 1.0, 1.1 & 1.2); both studies found the expected pattern of elevated R 1 values at the postero-medial aspect of Heschl's gyrus, reflecting putative auditory core. In visual cortical regions, Sereno et al. (2013) charted the increased R 1 values found within primary (e.g., V1) and non-primary (e.g., MT, V3a, V6) higher-visual areas, with respect to retinotopic functional borders (see also Abdollahi et al., 2014). Moreover, Sereno et al. showed overall higher R 1 in these visual regions compared to association cortex (angular gyrus), across intra-regional depths.
Here, we re-capitulate many of the findings above, and extend these results to further tissue-sensitive MPMs. Both R 1 and MT metrics in our present results manifested higher values within primary, primary adjacent, and higher visual cortical regions (e.g., Heschl's gyrus, probabilistic V1, probabilistic MT/V5, subcentral gyrus) compared to association areas (e.g., middle cingulate sulcus and gyrus, subparietal sulcus). In particular, that our inter-and intra-regional magnetization transfer (MT) results closely mirror our R 1 data lends strong support to the feasibility of mapping MT as a myelin proxy in the cortex. A potential benefit here is that MT is less affected by other properties of the tissue microstructure than R 1 -i.e., R 1 is partially influenced by R 2 * signal properties, which vary as a function of susceptibility effects due to paramagnetic molecule concentrations, myelin distribution and fibre orientationpoints we return to below (Callaghan et al., 2014b;Fukunaga et al., 2010;Stüber et al., 2014). MT therefore largely reflects the macromolecular content of the tissue microstructural environment, and can provide a proxy for the bound water fraction (Callaghan et al., 2014b). Moreover, previous post-mortem comparisons of MT ratio and T 1 -based MR metrics have shown high correlation between the two (r ¼ À0.79; Schmierer et al., 2004), suggesting that both converge well toward indexing myelin tissue. However, a difference in R 1 and MT in the current study was the reduced sensitivity of MT to myelin content within sulcal regions where cortex is thinnest, particularly across the calcarine fissure. As discussed in Results and above, the partial R 2 * influence on R 1 may have increased the contrast-to-noise for myelin and related processes within this sulcal region, where cortical thickness can fall to~1.5 mm, reducing the number of voxels that contribute to measured signal. The multiple signal contributions to R 1 therefore may enhance its sensitivity where measurements of tissue microstructure are required at a particularly fine-grained level, as in very thin and concave cortical regions.
A further consideration in the present results was the convergence between regions that showed high myelin content (i.e., elevated R 1 and MT values) and also high R 2 * values. The correspondence between R 1 , MT and R 2 * has been charted previously. In particular, tissue molecules such as iron influence R 2 * (Langkammer et al., 2010;Fukunaga et al., 2010) by causing local inhomogeneities in the B 0 field, in turn influencing the longitudinal relaxation rate (Callaghan et al., 2014b). Moreover, tissue regions high in iron often overlap areas of high myelin content, since oligodendrocyte cell bodies (whose cytoplasmic membranes extend around axons to form the myelin sheath) are known to express high concentrations of iron (Stüber et al., 2014;Bartzokis, 2004;Todorich et al., 2009). T 2 * (i.e., 1/R 2 *) has been shown to reflect patterns of myelin distribution, correlating highly with other myelin proxies in cortex (i.e., magnetization transfer ratio; Mangeat et al., 2015). T 2 * is also modulated as a function of fibre orientation, in line with measured fractional anisotropy (Cohen-Adad et al., 2012; see also Denk et al., 2011). Similarly, regions that show reduced MT in aging (reflecting presumptive demyelination) also tend to manifest reductions in R 2 * (Callaghan et al., 2014a), and demyelination in multiple sclerosis has been associated with focal increases in T 2 * as disease severity increases (Maneiro et al., 2015). In line with the inter-and intra-regional variation we observed in R 1 and MT, R 2 * tended to follow similar profiles of elevation or reduction within regions of respectively high or low R 1 and MT, mirroring its close relationship with R 1 and MT. Nevertheless, an important difference emerged in the location of high R 2 * foci, such that these were displaced into sulci adjacent to some of the regions of high R1 and MT (e.g., at Heschl's sulcus and gyrus, and at infero-temporal sulcus, adjacent to area MT; Supplemental Fig. 2). Further to the account of R 1 and MT differences (see 3.1, Group average MPM results), the displacement of R 2 * foci may reflect detection of iron-rich oligodendrocyte cell bodies within sulcal depths (e.g., Stüber et al., 2014). The exact mechanisms underlying this difference are unclear; however, a possible account is an interaction between the thinness of cortex over sulcal depths (as compared to thicker and more convex gyri), and the relative expression of macromolecular (i.e., myelin lipid) versus glial content as a result. The limited layer IV/V thickness within deep sulci constrains the expression of high macromolecular content, and in tandem, it is possible that oligodendrocytes may be more heavily expressed in sulcal regions proximal to heavily-myelinated gyri (see Stüber et al., 2014). In combination, these differences may lead to apparently elevated R 2 * in sulci adjacent to foci of high R 1 and MT.
Finally, we observed that intra-regional PD* parameters followed an expected pattern opposite that of R 1 , MT, and R 2 *, with highest values at cortical depths close to the pial surface. Agreeing with previous evidence of reduced myelin in superficial cortical layers (Annese et al., 2004;Walters et al., 2003;Leuze et al., 2014), the high PD* values likely reflect increased tissue water and lower myelination levels nearer the pial surface. The observed changes are in line with those that use [1 minus the ratio of PD and cerebrospinal fluid] as a marker for macromolecular tissue volume (Mezer et al., 2013).
An important consideration in interpreting MPM changes over cortical depth is the potential role of partial voluming effects. While it is possible that in thinnest regions of cortex, bounding depth fractions (i.e., 0.1 and 0.9) may be more susceptible to partial voluming, these effects would be expected to be constant over all MPMs, since the same surfaces are used per subject when mapping each parameter. Indeed, we found that cubic trends best explained the pattern of change in R 1 and MT over depth fractions, but a quadratic trend best explained the change in R 2 * over the same depths. Moreover, the apparent trends in the data are evident visually in Fig. 2, when restricting the range of depth fractions to 0.2-0.8a range we should expect to correspond with voxels inside the cortical sheet. If partial voluming effects drove the observed change over depth fractions, then cubic trends should be expected for each of the parameters, but this was not the case.

Development and cortical myelination
Here, we found evidence of a protracted course of myelin and myelinrelated process development within the cortex, across late adolescence and early adulthood. Our MPMs varied in the extent to which they revealed developmental effects (i.e., to hemisphere-wise FDR correction, q < 0.05): R 1 showed widespread increases with age across pre-frontal, frontal and parietal cortex; MT revealed a slightly lesser extent of age effects that were statistically significant with multiple comparison correction, largely concentrated in pre-central regions bilaterally. Age effects on R 2 * were not robust in vertex-wise analyses, but did emerge in multi-level models that included data from each depth fraction per ROI.
Of particular relevance to our present results, recent studies using T 1w /T 2w image contrast ratios have identified extended periods of change within association cortex. Shafee et al. (2015) and Grydeland et al. (2013) found significant increases in T 1w /T 2w ratio with age over much of the frontal and parietal lobes, in young adults (18-35 years) and across the lifespan (8-83 years), respectively. Notably, the linear trend identified in the present study (and by Shafee et al., 2015) was also found by Grydeland et al. (2013) when considering the younger (8-20 year old) tail of their age distribution (cf. quadratic trends in T 1w /T 2w ratio across their full age range). This appears to reflect phases of increasing myelination up to early adulthood, followed by periods of relative stability and eventual decline of cortical myelin after 50-60 years of age (Grydeland et al., 2013; see also Miller et al., 2012).
As outlined above, an important advance made by our age findings is our use of MPMs that index a range of myelin-related tissue processes. MPMs enabled us to measure tissue-based parameters reflecting water mobility, macromolecule concentration, and paramagnetic ions (R 1 ), together with either predominantly macromolecular (MT), or paramagnetic ion (i.e., iron) (R 2 *) concentration (Callaghan et al., 2014b). Existing studies that have used T 1w /T 2w ratio methods cannot resolve for tissue-sensitive parameters, since both T 1w and T 2w contrast are determined by a variety of microstructural (see Glasser et al., 2014;Vidal--Piñeiro et al., 2016), and vasodilatory properties (i.e., O 2 /CO 2 concentration; Tardif et al., 2017). Our finding of regions that manifested overlapping age effects for R 1 and MT suggests that myelination follows a protracted developmental course. Interestingly, we did not observe robust vertex-wise age effects on R 2 *, a parameter sensitive to tissue iron concentration, myelination and fibre orientation. Whereas previous studies have shown focal R 2 * decreases along fibre bundles (cf. sub-cortex) in line with aging in older samples (Callaghan et al., 2014a), here, our younger cohort showed most robust increases in parameters largely reflective of cortical myelin content (although we note multi-level analyses of ROI data did reveal age effects for iron and fibre-orientation sensitive R 2 *). One speculative account is that the developmental effects we observed involved changes to myelin sheath thickness (i.e., g-ratio; Dean et al., 2016) rather than large-scale changes in the numbers or density of iron-rich oligodendrocyte cell bodies. However, optogenetic evidence in mice has supported effects of behaviorally-relevant neural activity in promoting increases in both myelin sheath thickness (i.e., g-ratio decreases) and oligodendrogenesis (Gibson et al., 2014). Taken together, it is likely that our current effects of age reflect some combination of these processes, although the precise mechanisms remain unclear. Future studies in which in-vivo measurements of g-ratio (e.g., Mohammadi et al., 2015) are probed across our age range in addition to each multi-parameter map may enable us to shed further light on the mechanisms at play in myelin sheath development.
A strength of the present approach is the use of a well-validated MPM protocol incorporating B 0 and B 1 RF transmit field mapping. Here, field maps form an integral part of the MPM protocol, such that local flip angles can be resolved and used in the estimation of MPMs, based on a variable flip angle procedure (see Lutti et al., 2010Lutti et al., , 2012Helms et al., 2008b;Weiskopf et al., 2011). T 1w /T 2w ratio methods are subject to B 1 RF transmit field inhomogeneities during acquisition, which can bias local flip angles and the resulting signal intensity and contrast, increasing measurement error (Lutti et al., 2010Helms et al., 2008b;Weiskopf et al., 2011; but see also Glasser et al., 2014). The precision and reproducibility of the MPM and field mapping protocols have been documented previously Lutti et al., 2010Lutti et al., , 2012. Importantly, although we observed differences in MPM parameter values between the protocols that differed across cohorts in our present study, those differences reflected a constant offset, the source of which was isolated and which we controlled for in statistical models. A potential limitation of our present R 2 * measurements arises from the fact that we did not correct subject-wise R 2 * maps for local fibre orientation relative to the B 0 field. Previous investigations have shown that the difference in angle between vertex-wise surface normal projections and the B 0 field can account for some of the measured T 2 */R 2 * in heavily myelinated regions (Cohen-Adad et al., 2012;see Denk et al., 2011). Nevertheless, the extent of variation captured by B 0 dependence as previously reported would correspond with <5% of R 2 * based on the midpoint presented in our Fig. 1. Thus, we suggest that in large part, our present reported values are representative of much of cortex, with the caveat that R 2 * values from some heavily myelinated regions (e.g., BA4) may be subtly lower than the values reported here if corrected for fibre orientation. The use of sufficiently high-resolution diffusion data (not acquired as part of this study) may help to validate any B 0 dependence with respect to intra-cortical fibre orientation.
Nevertheless, in light of the methodological advances in MPM methods, our present findings agree well with accounts of cortical development and myelination based on T 1w /T 2w ratio methods. Moreover, that we were able to identify robust linear developmental effects on cortical myelination/myelin processes and with a much smaller cohort than many studies (Shafee et al., 2015;Grydeland et al., 2013) suggests that reliable quantitative mapping protocols can play a highly informative role in charting cortical development.

Conclusions
Using a well-validated multi-parameter mapping protocol, we showed that both inter-and intra-regional cortical myelin and related processes can be quantified across much of the cortical sheet. We further demonstrated the utility of delineating profiles of cortical myelin processes across cortical depths (using R 1 , MT, and R 2 *) in tandem with mapping of tissue water sensitive parameters (PD*). Moreover, exploring effects of development cross-sectionally, we found that cortical myelin and myelin processes increased at a rate of 0.5 %-1.0% (R 1 and MT) per year over frontal and parietal regions, across late adolescence and early adulthood. These results shed further light on ontogenetic factors that may shape large-scale cortical organization and inform broader accounts of lifespan cortical development, as well as helping to characterize the healthy aging of the human brain, which may provide a useful clinical benchmark for studying de-myelination in aging and disease.

Funding
The research leading to these results received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement n 616905, and via EC FP7 grant n MC-ITN-264301 (TRACKDEV) to DC. The Wellcome Trust Centre for Neuroimaging is supported by core funding from the Wellcome Trust 0915/Z/10/Z. The work was also supported by a Wellcome Trust grant 100227 (M.A., G.R.).

Conflicts of interest
The authors declare no competing financial interests.