Depth-dependent intracortical myelin organization in the living human brain determined by in vivo ultra-high field magnetic resonance imaging

BACKGROUND
Intracortical myelin is a key determinant of neuronal synchrony and plasticity that underpin optimal brain function. Magnetic resonance imaging (MRI) facilitates the examination of intracortical myelin but presents with methodological challenges. Here we describe a whole-brain approach for the in vivo investigation of intracortical myelin in the human brain using ultra-high field MRI.


METHODS
Twenty-five healthy adults were imaged in a 7 Tesla MRI scanner using diffusion-weighted imaging and a T1-weighted sequence optimized for intracortical myelin contrast. Using an automated pipeline, T1 values were extracted at 20 depth-levels from each of 148 cortical regions. In each cortical region, T1 values were used to infer myelin concentration and to construct a non-linearity index as a measure the spatial distribution of myelin across the cortical ribbon. The relationship of myelin concentration and the non-linearity index with other neuroanatomical properties were investigated. Five patients with multiple sclerosis were also assessed using the same protocol as positive controls.


RESULTS
Intracortical T1 values decreased between the outer brain surface and the gray-white matter boundary following a slope that showed a slight leveling between 50% and 75% of cortical depth. Higher-order regions in the prefrontal, cingulate and insular cortices, displayed higher non-linearity indices than sensorimotor regions. Across all regions, there was a positive association between T1 values and non-linearity indices (P < 10-5). Both T1 values (P < 10-5) and non-linearity indices (P < 10-15) were associated with cortical thickness. Higher myelin concentration but only in the deepest cortical levels was associated with increased subcortical fractional anisotropy (P = 0.05).


CONCLUSIONS
We demonstrate the usefulness of an automatic, whole-brain method to perform depth-dependent examination of intracortical myelin organization. The extracted metrics, T1 values and the non-linearity index, have characteristic patterns across cortical regions, and are associated with thickness and underlying white matter microstructure.


Introduction
Nearly half the human brain comprises myelinated axons that connect brain regions with each other. Myelin is formed by oligodendrocytes that sheath and insulate neuronal axons (Nave and Werner, 2014). Myelination supports the speed of axonal conduction which determines the time of arrival of signals from multiple axons and hence neuronal spiking (Nave and Werner, 2014). Myelin also provides trophic axonal support which is particularly relevant for longer axons that extend far beyond the neuronal soma (Bercury and Macklin, 2015). Myelination is therefore essential for the efficient inter-neuronal communication that underlies mental processing and experience-induced plasticity (Fields et al., 2014;Fields, 2015;Hunt et al., 2016;Glasser et al., 2014).
Although myelination is a prominent feature of deep white matter, oligodendrocytes are also abundant within the cortical ribbon. Within the cortex, myelin density is higher in deeper compared to superficial layers and in sensorimotor compared to association regions (Nave and Werner, 2014). Intracortical myelin has a prolonged period of maturation extending into late adolescence and adulthood (Carey et al., 2017;Shafee et al., 2015;Whitaker et al., 2016) which is associated with age-related changes in gray matter volume (Mills et al., 2016) and in functional connectivity (Huntenburg et al., 2017). Intracortical myelin is sensitive to experience-dependent neuronal activity throughout the lifespan thus actively contributing to brain plasticity and remodeling (Fields, 2015;Purger et al., 2016). These features suggest that myeloarchitecture is relevant to the fine-tuning of cortical circuits, a notion that is supported by the association between intracortical myelin and individual variability in cognitive function (Grydeland et al., 2013).
The majority of studies on cortical myelination to date have been performed at the standard field strength of 3 Tesla with typical spatial resolutions of 1 mm. Since the thickness of the human cortex ranges from approximately 2 mm-4 mm, these standard MR techniques are suboptimal for the characterization of intracortical myelin concentration. Ultrahigh field imaging at 7 Tesla now allows for the sub-millimeter resolution in MR sequences (Fracasso et al., 2016;Waehnert et al., 2016;Marques et al., 2010Marques et al., , 2017Sanchez-Panchuelo et al., 2010;Tardif et al., 2015). Initial studies suggest that depth-dependent intracortical myelin properties in the sensorimotor  and the visual (Fracasso et al., 2016) cortices obtained using ultra-high field MRI are consistent with post mortem histological data. MRI-based laminar T 1 maps align with cortical folding patterns ) and with regional functional specialization (Sanchez-Panchuelo et al., 2010;Kuehn et al., 2017), further supporting a close relationship between plasticity and laminar myelin.
The ability to approximate the distribution of myelin across lamina by sub-sampling the cortex at different depth-levels is an important new development for two reasons. Firstly, myelin is highly unequally distributed across cortical lamina. Hence, measuring myelin at low spatial resolutions yields metrics that are influenced by partial volume effects across lamina, and are therefore less sensitive to inter-regional or inter-individual variation than depth-specific estimates. Secondly, cortical myelination across development (Konner, 1991) and the myelination response to inflammation (Gardner et al., 2013) are lamina-specific processes. Thus, myelin concentration in different cortical layers can play distinct roles in neuronal circuit functioning, neurodevelopment and pathology. In psychiatric disorders such as schizophrenia, the importance of layer-specific cortical investigations has also been emphasized in the context of histological studies (H€ oistad et al., 2013;Rajkowska et al., 2002).
In order to study lamina-specific myelin pathology in histology, each layer tends to be investigated separately for clinical associations with cell density or myelin stains. However, explorative whole-brain MRI studies would benefit from a single marker for cross-laminar distribution of myelin, as increased multiple testing results in low statistical power. In the present study, we present an optimized whole-brain MR sequence and automated pipeline for the investigation of the myeloarchitecture of the brain in healthy individuals using ultra-high field imaging. In addition to analyzing T 1 values at each cortical depth level, we explored myelin profile nonlinearity as a single index of myelin organization. To aid interpretation of this new metric, we provide information on how this nonlinearity index relates to standard MRI measures of cortical thickness and sub-cortical white matter fractional anisotropy (FA) measured using diffusion weighted imaging.

Sample
The sample consisted of 25 healthy volunteers (14 male and 11 female) with a mean age of 29 years (sd ¼ 5.88) recruited through adverts in the local community. They were interviewed to exclude any past or current medical disorder or head trauma and any psychiatric disorder based on the Mini-International Neuropsychiatric Interview (M.I.N.I.) (Sheehan et al., 1998). We also recruited 5 patients from the local specialist clinic diagnosed with multiple sclerosis (MS) (2 male and 3 female; age range: 32-55 years), as "positive controls" since MS represents the prototypical demyelinating condition (details in Supplement). The study was approved by the institutional review board of the Icahn School of Medicine at Mount Sinai (ISMMS). All study participants provided written informed consent.

MRI acquisition
All participants were scanned at ISMMS using a 7T MR scanner (Magnetom, Siemens Healthcare) with a 32-channel with a Nova head coil (Nova Medical, Wilmington MA). Whole brain T 1 -weighted images were acquired using ultra-high resolution MP2RAGE sequence with the following parameters: 0.5 mm isotropic resolution, repetition time (TR) ¼ 5000 ms, echo time (TE) ¼ 5.75 ms, inversion time (TI) TI 1 / TI 2 ¼ 900 ms/2780 ms, 224 axial slices with slab thickness 11.5 cm, fieldof-view ¼ 224.5 Â 203 Â 112 mm 3 , and slab selective excitation and flow suppression (Marques et al., 2010). Total scan time was 25 min. Images at the two inversion times (TI 1 /TI 2 ) were used to calculate bias-field corrected quantitative T 1 maps and T 1 weighted images.
Diffusion-weighted images were also acquired in the same session with a single shot spin echo EPI sequence with monopolar diffusion encoding Vu et al., 2015). Parameters included: 2 mm isotropic resolution, 108 directions/b-values from 0 to 2000 s/mm 2 , TR/TE ¼ 4000/62 ms, whole brain coverage, multi-band factor of 2 , in-plane GRAPPA acceleration R ¼ 3. Paired acquisitions with reversed phase encoding in the AP/PA direction were acquired at 7.5 min per scan for a total scan time for dMRI of 15 min. The gradient table for the scan was identical to that presented in Chiang and colleagues (Chiang et al., 2014) with 99 diffusion-encoding directions isotropically distributed from b ¼ 200 s/mm 2 to b ¼ 2000 s/mm 2 with 9 b ¼ 0 s/mm 2 images interleaved. De-identified data may be made available to researchers in other institutions, upon request to the corresponding author, under a legally binding data use agreement that specifies of the planned data use.

High resolution intracortical myelin profiling using T 1 maps
Image preprocessing followed the procedures developed by Dinse and colleagues . The T 1 maps, T 1 -weighted images and TI 2 images were skull-stripped, background-masked and aligned (using rigid 6 degrees-of-freedom registration, normalized mutual information) to a 0.4 mm MNI template using the CBS Tools (https://www.nitrc.org/ projects/cbs-tools). Next, images were corrected for field inhomogeneities and histogram-matched to a template image from the sample of  using 3D Slicer (https://www.slicer.org/). Using the Topology-preserving, Anatomy-Driven Segmentation (TOADS) (Bazin and Pham, 2007) and Multi-object Geometric Deformable Model (MGDM)  segmentation algorithms as integrated in the CBS Tools, images were filtered with arteries and dura segmentations, corrected for partial volume effects and were then segmented into gray matter, white matter and cerebrospinal fluid (CSF). The cortex was extracted using the CRUISE algorithm (Landman et al., 2013), which is robust to white matter lesions (Shiee et al., 2014), and divided into 20 cortical depth-levels using a volume-preserving approach as implemented in the Volumetric Layering module of CBS Tools (Waehnert et al., 2014). This level-set approach generates one 2-dimensional surface per layer containing T 1 values at each surface location. Next, the Profile Sampling module of CBS Tools was used to transform the 2-dimensional level-sets to 3-dimensional representations of T 1 values. This transformation represents each vertex T-values of each level-set surface as a column of identical T 1 values perpendicular to the surface . The final result of this pipeline is 20 vol of cortical profiles, one for each depth-level, for each of 15 healthy and the 5 MS patients. The whole process was completed within each individual's native space. A set of exemplary results from healthy individuals and patients with MS are shown in Supplemental Figures 1 and 2.
Next, the cortex was segmented and parcellated in native space using Freesurfer (https://surfer.nmr.mgh.harvard.edu/). This parcellation was based on the pre-processed T 1 -weighted images alone (independent of the cortex extraction described above) and was used to obtain 148 cortical regions-of-interest (ROIs) in native-space based on the Destrieux atlas. We used the Destrieux atlas, which separates gyri from sulci (Destrieux et al., 2010) because cortical myelin profile and content varies between gyri and sulci (Waehnert et al., 2014;Sereno et al., 2013). Prior to further analyses, we implement a strict quality control protocol whereby whole brain 3D segmentations were visually inspected and we further checked the Freesurfer cortical thickness estimates for regions around long/short gyri and the claustrum where segmentation is often challenging. Mean thickness values, for left and right respectively, were 3.44 mm and 3.37 mm for the short gyrus, 2.92 mm and 2.94 mm for the anterior insula, 2.59 mm and 2.50 mm for the inferior insula, and 2.54 mm and 2.74 mm for the superior insula. Overall mean cortical thickness according to FreeSurfer was 2.47 (sd ¼ 0.32).
In each participant and for each cortical ROI, we estimated (a) the absolute T 1 value as an index of myelin concentration at each depth-level, (b) the non-linearity index, a novel measure reflecting the depthdependent distribution of T 1 values along the cortical depths from the outer CSF boundary to the inner white-matter boundary, (c) distinct positive and negative deviations from linearity, (d) cortical thickness as provided in the Freesurfer output. The non-linearity index is the deviation from linearity from the curve of T 1 value distribution along the cortical ribbon, which was calculated as the squared difference between the observed T 1 values and the linear regression line approximated by the curve (Fig. 1). The positive and negative deviations from linearity represent respectively the area above and underneath the regression line in Fig. 1. Fig. 1 also contains two illustrative examples from occipital and prefrontal regions, respectively with low and high profile nonlinearity. We chose to derive cortical thickness from Freesurfer on the basis of the T 1 -weighted maps alone, independent of the myelin profiling pipeline, to minimize confounding of cortical thickness estimation with the T 1 values and the cortex extraction. The scripts we used to run the pipeline is available upon request, here.
For the MS patients, we simply transformed their T 1 values to Z-scores based on the T 1 distributions in the healthy individuals, and inspected the degree of deviation from the normative values as a qualitative validation of the interpretation of the T 1 profiles.

Quality assessment of T 1 profiling in the healthy participants
Several steps were undertaken to ensure the data quality of the final T 1 profiles. First, three participants were excluded after visual inspection of the raw volumes and the Freesurfer cortical parcellations. Second, across all participants we excluded individual regions with outlier values based on the deviation from the median T 1 profile of each region. To this end, we calculated shape of the median profile of each region which closely a cubic spline (all R 2 > 0.98). Individual regions were excluded case-wise (1) their fit (R 2 ) to the median spline was <0.70; or (2) their absolute deviation from the median spline was at the 5% extreme end of the sample. These criteria allowed stricter control with regards to the shape of the cortical myelin profiles and less strict control on the height of the curve. This is because we considered significant deviations from the typical regional myelin profile shape as indicators of poor data quality. A deviation in myelin content at a specific depth-level may, on the other hand, indicate subtle but real differences in cortical organization. Following the application of these quality control criteria, 34 of the 148 cortical regions were excluded completely from the analysis because many participants had deviant profiles (Supplemental Table 1). Figures of all myelin profiles are available in Supplemental Information-2, with profiles excluded after quality control depicted in red. Visual inspection of these figures suggests that the objective quality control criteria was adequate at consistently identifying profiles that deviated a lot from typical patterns. In accordance with previous reports (Glasser and Van Essen, 2011;Haast et al., 2016), artifacts due to residual inhomogeneity were observed in the medial orbitofrontal and inferior temporal lobes. We also observed artifacts in the most superior medial fronto-parietal regions due to RF inhomogeneity and reduced field of view. As mentioned, these regions were removed from the analyses following our quality control procedures. The data for the MS patients were subject to quality control as described above but we did not exclude patients on the grounds of deviation the myelin profiles because such deviations from "typical" profiles are likely in this population. Generally, it is (as always) possible that some subtle artifacts or segmentation biases, not visible by the naked eye or not captured by our quality control algorithm, could have persisted. However, to really alter our conclusions, they would have to be non-random in relation to cortical depth and thickness, and this we consider unlikely.

Extraction of sub-cortical white matter FA values
Diffusion-weighted MRI data were pre-processed in accordance to pipelines developed by the Human Connectome Project (HCP)  modified in-house to accommodate the imaging protocol and the B0-eddy currents of our scanner (O'Halloran et al., 2015). Multi-band dMRI images were reconstructed using SENSE . Preprocessing consisted of slice-wise eddy current correction of the first 2 acquired slice groups, fitting of the field map from oppositely encoded non-diffusion-weighted (b ¼ 0) images with "topup" in FSL, correction of distortions due to eddy currents using FSL, and correction of gradient non-linearities using "gradunwarp", an HCP routine implemented in python. The diffusion-weighted volumes were visually inspected for quality assurance. The data were then skull-stripped and diffusion tensors were calculated respectively using "bet" and "dtifit" in FSL. To extract white matter FA-values underneath each cortical region, the Destrieux cortical regions were dilated with 1 mm, binarized, and multiplied by the white matter segmentation. For each subject, the individual T 1 -weighted image was linearly transformed to each subject's b 0 image using "FLIRT2 in FSL with 6 degrees of freedom and used to transform the Freesurfer ROIs to the FA map. The subcortical regional masks were then overlaid on the FA maps to extract the mean FA values underneath each region (Supplemental Figure 3).

Statistical analyses
Statistical analyses were implemented in R. To help interpret the nonlinearity index, we first investigated the contribution of each depth-level to the non-linearity index. For this, we used stepwise mixed linear models with non-linearity as dependent variable, cortical region as fixed factor, subject as random factor (with random intercept), and all T 1 values at each depth level as separate fixed regressors. Next, analyses of T 1 values (indicative of myelin concentration) and analyses of the non-linearity index were conducted separately using the same mixed linear models with cortical thickness as continuous predictor. We also calculated Spearman's rank correlation coefficients between cortical thickness and T 1 values and non-linearity indices within each region across all individuals as well as within each individual. The aim of these correlational analyses was to confirm that the pattern identified using the mixed regression models was also observed at the level of univariate correlations and not to make inferences about statistical significance. The same methods (mixed-models and correlations) were used to investigate associations between subcortical FA values with T 1 and non-linearity indices. To investigate whether positive deviations or negative deviations were mainly driving the observed associations with nonlinearity, significant associations of the non-linearity index were also tested for positive and negative deviations from linearity specifically, using the same mixed linear model approach.
We then tested whether regions formed separate clusters with respect to their average non-linearity using Gaussian mixture modeling clustering as implemented in the mclust package in R (Scrucca et al., 2016). This algorithm uses an expectation-maximization algorithm to decide upon the optimal clustering solution and also provides the Bayesian information criterion (BIC) to corroborate the choice.

Properties of intracortical myelin profiles in healthy individuals
As expected T 1 values were (a) lower near the gray-white matter boundary and higher near the CSF surface with a slight leveling between 50% and 75% of cortical depth and (b) lower around the central sulcus and in the visual and auditory cortices. Fig. 2A shows the median T 1 values at three cortical depth-levels. At each level, the expected increase in myelin around the central sulcus, the visual cortices and the auditory cortices are clearly visible. Fig. 2B shows the distribution of the nonlinearity index across the cortical ribbon. Although cortical regions differed in their non-linearity index these differences did not form discreet clusters (Supplemental Table 2). Stepwise mixed model regression of nonlinearity across all regions, with T 1 values of all depth levels as predictors revealed that, in general, T 1 values at 5%, 10%, 80% and 100% were the strongest negative predictors of non-linearity, whereas T 1 values at 40%, 55% and 65% were the strongest positive predictors or profile nonlinearity (all P < 0.0025). This indicates that nonlinearity across regions was driven mostly by a T 1 "plateau" of lower myelin in the middle layers in combination with higher myelin estimates at ends of the curve. The nonlinearity indices of outliers were differently influenced by T 1 values at each depth-level (interactions with depth-level: 0.01 < P < 10 16 ), indicating that nonlinearity profiles of these outliers do not always represent the same characteristics as for the typical profiles.  Figure 4 show individual MS patients' deviation from controls in terms of overall and region-specific intracortical T 1 values. There was a high degree of inter-individual heterogeneity. One patient very clearly had widespread and pronounced (Z > 4) reductions in estimated myelin concentrations. At the other extreme, one patient appeared to have increased myelin concentration in several, different cortical regions. Taken together, results were most consistent at deeper cortical levels around 75% of depth, and around the central sulcus and in the cingulate cortex. In both of these regions 3 of the five patients showed myelin levels that were 2-5 standard deviations lower than those in controls.

Association between T 1 values and cortical thickness
Across all regions (modeling region as a fixed factor and subject as a random factor), there was a positive association (4.23 < t < 14.69; P values between 10 À4 and <10 À13 ) between cortical thickness and T 1 values at deeper (30%-100%) but not outer (0%-25%) cortical depthlevels (Supplemental Table 3). This pattern was also apparent in the univariate correlations across and within individuals. The correlation between T 1 values and cortical thickness was positive across individuals at 75% of cortical depth but not at 25% (Supplemental Table 4). Within individuals, T 1 values across regions also correlated consistently with cortical thickness in the deeper (ρ range: 0.32 to 0.62) rather than outer cortical depth-levels (ρ range: À0.10 to 0.38) (Supplemental Table 5). There were no noteworthy effects of age, age 2 but there were region-bysex interactions on T 1 at 25% (P ¼ 0.003), 50% (P ¼ 5.83*10 À5 ), and 75% of cortical depth-level (P ¼ 1.88*10 À7 ); women tended to have higher T 1 values compared to men (Supplemental Table 6).

Association between non-linearity indices and cortical thickness
Across all regions (modeling region as a fixed factor and subject as a random factor), there was a positive relationship between non-linearity indices and cortical thickness (t ¼ 8.27; P < 10 À15 ). This was a consistent observation across and within individuals (Supplemental Tables 7  and 8). The association between cortical thickness and profile nonlinearity was equally driven by both positive (t ¼ 6.30; P < 10 À9 ) and negative (t ¼ 6.33; P P < 10 À9 ) deviations from linearity. Similar analyses also found a positive relationship between non-linearity index and average T 1 values (t ¼ 4.54; P < 10 À5 ). The prefrontal and cingulate cortices and the insula, displayed markedly more profile nonlinearity than sensorimotor regions such as occipital regions, the superior temporal gyri, and the postcentral gyri (Fig. 2B). When we included the outliers in the analysis, the overall association of nonlinearity with cortical thickness in the general mixed model remained significant, although somewhat attenuated (all P < 10^-5). There were no noteworthy effects of age, age 2 or sex on non-linearity indices in any region (all P > 0.39).

Association between subcortical FA and intracortical myelin measures
There was an interaction between depth-level and FA on T 1 value (F ¼ 80.41, P < 2.2*10 À16 ). At 95% of cortical depth every individual displayed a negative correlation between FA and T 1 value, ranging from À0.13 (P ¼ 0.18) to À0.46(P < 10^-7) (Fig. 4), indicating that regions with higher myelin concentration at the gray matterwhite matter boundary also have higher FA below this boundary. Correlations at multiple cortical depth-levels, within regions, and within individuals, are available in Supplementary Tables 10-11. There was a positive correlation between subcortical FA and non-linearity indices across regions for all individuals (0.06 < ρ < 0.54; 10 À9 < P < 0.50, Supplemental Tables 12-13). For both positive and negative deviations from linearity correlations with FA were between 0.10 (P ¼ 0.29) and 0.56 (P < 10-9). When we included the outliers, the individual correlations with FA remained consistently positive (0.10 < rho < 0.54).

Discussion
We used quantitative ultra-high field MR imaging to examine intracortical myelin organization in the living human brain in relation to cortical thickness and subcortical FA values. In addition to estimating myelin concentration, we computed the non-linearity index, a novel metric that assesses the spatial distribution of myelin along the cortical ribbon.
We found that estimates of myelin concentration were higher in sensory and motor cortices and lower in higher-order association areas ( Fig. 2A). Our results confirm previous findings in healthy individuals obtained using the same acquisition and processing methods  as well other methods such as T 1 /T 2 ratios (Glasser et al., , 2016Haast et al., 2016;Marques et al., 2017;Van Essen and Glasser, 2014), and T 2 * (Haast et al., 2016;Marques et al., 2017). Qualitative comparison with 5 patients with MS, used as "positive controls", showed that the current method identified cortical myelin concentration alterations in this prototypical demyelinating disease.
Cortical thickness and estimated myelin concentration (inferred from T 1 values) were positively associated. This association was present across regions but also across individuals. This interregional variability in cortical thickness and myelin concentration is also thought to underlie inter-individual variability (Grydeland et al., 2013) although methodological explanations relating to cortical segmentation cannot be fully ruled out. Failure to include highly myelinated deep cortical layers would result in lower myelin concentration estimates while erroneous inclusion of subcortical white matter would lead to higher estimates. However, our results are less likely to be confounded by segmentation errors given that the cortical thickness measures were derived from a different segmentation than the myelin estimates.
The results reported here suggest that subcortical FA values are associated with higher myelin concentrations at the deeper cortical levels (Fig. 4). FA reflects a number of white matter properties including, but not limited to, myelin. In general, FA values estimated near the graywhite matter boundary are more susceptible to artifacts in the modeling of white matter fibers close to the cortex. However, regions where highly myelinated fibers enter or leave the cortex are expected to have both higher FA and lower T 1 values. The consistency of this association across subjects, and its lack in the outer cortical layers, supports a neurophysiological interpretation. Moreover, confounding effects due to cortical segmentation errors are unlikely, because over-estimation of the cortical depth is less likely for high-FA, stereotypical white matter. Thus, the association between sub-cortical FA and quantitative T 1 values is in line with neurophysiological properties and represents a cross-modal confirmation of the T 1 profiling approach used here.
A key novelty of this paper is the introduction of the non-linearity index as a measure of intracortical myelin organization (Fig. 1). The availability of an automatically generated, single metric for myelin organization obtained from in vivo MRI benefits future, larger-scale, datadriven studies in terms of statistical power and required resources. In addition, in studies of pathology, the nonlinearity index will be more sensitive to deviations from the norm in the presence of etiological heterogeneity in the type of myelin disorganization across patients, as may be the case in major psychiatric disorders. Variability in the non-linearity index was primarily driven by decreased myelin in the middle cortical layers, in combination with increased myelin at both outer and inner depth-levels. Associations of total nonlinearity with other anatomical metrics was driven by equal contributions of both decreased myelin in middle layers and increased myelin at both inner white matter and outer CSF boundaries. The interregional correlations of T 1 values and profile nonlinearity with cortical thickness and FA likely reflect general organizational principles of the brain linking regional morphology and connectivity (Barbas, 2015;Beul et al., 2015;Huntenburg et al., 2017;Scholtens et al., 2014). For example, unimodal sensory cortical are relatively thin, have more myelinated neurons, fewer and less elaborate dendrites, and sparser functional connections compared to higher-order areas (Beul et al., 2015;Huntenburg et al., 2017;Scholtens et al., 2014). Higher-order regions, on the other hand, contain less myelin, have more abundant connections with the rest of the brain (Beul et al., 2015;Huntenburg et al., 2017;Scholtens et al., 2014), and more complex Fig. 4. Fractional anisotropy (FA) in the segment just below the cortex was negatively associated with T 1 values in deep cortical layers. Each line represents one subject, showing the consistency of these relationships across individuals. Each line is a linear fit of one subject's T 1 values at 95% depth regressed on subcortical FA values across the brain. dendritic arborization that takes up more unmyelinated space within the cortex (Glasser and Van Essen, 2011). These organizational principles are also reflected in distinct evolutionary origins and developmental trajectories, where ontological and phylogenic early regions serve more primary functions and are less plastic throughout life, compared to late evolved and late developing cortical regions. Raw T 1 values, cortical profiles and nonlinearity indices appear to capture partly overlapping and partly unique aspects of these general organizational principles of the brain.
Taken together, our results suggest that the non-linearity index effectively captures inter-regional and inter-individual variation in laminar myelin organization that is not represented by absolute T 1 values or other metrics alone. However, an understanding of the physiological properties underlying the non-linearity index is beyond the resolution of neuroimaging. Evidence from post-mortem human studies suggests that the higher non-linearity index observed in prefrontal regions may be due to lower myelin content in deep layers (Nave and Werner, 2014). On the other hand, investigations of the mouse cortex have shown that in the superficial cortical layers 2 and 3 there exists a surprising amount of a distinct type of myelin surrounding inhibitory neurons (Micheva et al., 2016), raising the possibility that an increase in this myelin could influence nonlinearity. The non-linearity index may also be influenced by cortical thickness per se as thicker cortices allow for improved sampling of the T 1 values across the cortical ribbon leading to smaller partial volume effects. Post-mortem evidence in humans and animal studies are required to better understand the physiological underpinnings of the non-linearity index. Given the importance of intracortical myelin organization for human behavior and neuropsychiatric disorders (Haroutunian et al., 2014), the non-linearity index may be an enticing new imaging phenotype to investigate in clinical populations and across the lifespan. A notable advantage of the present work is that is extends previous manual procedures of cortical myelin profiling to become automated and data-driven, and therefore applicable to large datasets and exploratory research in bran disorders with unknown etiology.
Several limitations of this study ought to be clarified. Firstly, the signal-to-noise ratio was relatively low, but still acceptable, due to the high resolution of our quantitative whole-brain T 1 maps. The inhomogeneous transmit field at 7T caused considerable signal dropout in the most inferior regions and frequent clipping of the superior medial regions. These cortical regions were therefore removed from our analysis (Supplemental Table 1). This tradeoff allowed us to achieve maximum brain coverage at 0.5 mm isotropic resolution in a reasonable scan time ensuring that the data acquisition would be easily tolerated by participants. Issues of tolerability and cost-containment are particularly important for the application of this method to larger samples and to clinical populations. We choose to use cortical parcellation, rather than vertex-wise analyses, to avoid normalization to a common template, which requires reshaping individual cortical surfaces to a template and can obscure inter-individual differences in folding. Moreover, we used the Destrieux cortical parcellation scheme because it was developed according to precise anatomical rules, taking into account cytoarchitectural boundaries and curvature landmarks, and it produces anatomical labels that are closely matched to expert manual tracings (Destrieux et al., 2010). Therefore, whilst we might lose some anatomical resolution, we argue that our approach achieves greater power, minimal confounding of inter-individual variability in curvature and improved interpretability, compared to template-normalized vertex-wise analyses. We note the obvious limitation that T 1 values do not directly reflect myelin concentrations, but are a close approximation. T 1 values are also influenced to a lesser extent by iron concentrations although this is unlikely to be a major factor in a sample of young healthy individuals. Despite these limitations, we obtained quantitative measurements that passed quality control, and are consistent with previous in vivo and ex vivo reports from independent research groups using variable methods. Intracortical myelin quantification and high-field MRI are very fast developing areas of research, and the present findings contribute to this work in progress.
In conclusion, we demonstrate the use of an automatic, whole-brain method to perform data-driven examination of intracortical myelin organization using quantitative MRI in vivo. The extracted metrics, T 1 values and the non-linearity index, have characteristic patterns across cortical regions, and are associated with thickness and underlying white matter microstructure. This approach enables large-scale and longitudinal investigations of intracortical myelin throughout development and aging and in clinical populations.

Funding
This work was supported in part through computational resources and staff expertise provided by Scientific Computing at the Icahn School of Medicine at Mount Sinai, USA. Dr. Frangou received support from the National Institutes of Health, USA (R01MH113619). Dr. Sprooten is supported by a Hypatia Tenure Track grant from the Radboud University Medical Center, Netherlands and a NARSAD Young Investigator Grant from the Brain and Behavior Foundation, USA (ID: 25034).