Mild Traumatic Brain Injury Results in Significant and Lasting Cortical Demyelination

Despite contributing to neurocognitive deficits, intracortical demyelination after traumatic brain injury (TBI) is understudied. This study uses magnetic resonance imaging (MRI) to map intracortical myelin and its change in healthy controls and after mild TBI (mTBI). Acute mTBI involves reductions in relative myelin content primarily in lateral occipital regions. Demyelination mapped ~6 months post-injury is significantly more severe than that observed in typical aging (p < 0.05), with temporal, cingulate, and insular regions losing more myelin (30%, 20%, and 16%, respectively) than most other areas, although occipital regions experience 22% less demyelination. Thus, occipital regions may be more susceptible to primary injury, whereas temporal, cingulate and insular regions may be more susceptible to later manifestations of injury sequelae. The spatial profiles of aging- and mTBI-related chronic demyelination overlap substantially; exceptions include primary motor and somatosensory cortices, where myelin is relatively spared post-mTBI. These features resemble those of white matter demyelination and cortical thinning during Alzheimer's disease, whose risk increases after mTBI.

Whereas both acute and chronic post-traumatic demyelination of the cerebral cortex have been studied extensively in animals (9)(10)(11), hardly any human studies have quantified the spatial profiles and extent of these phenomena. This is due partly to the challenges of measuring human cortical myelin content in vivo. The gold standard for measuring human cortical myelin content is postmortem histopathological examination. However, as Glasser & Van Essen (12) showed, the ratio R of T 1 -to T 2 -weighted magnetic resonance images (MRIs) provides an accurate measure of brain myelin content because T 1 -weighted MRI intensity and the reciprocal of T 2 -weighted MRI intensity are both proportional to myelin content. Computing R offers improved myelin estimates because this measure implicitly adjusts for nuisance effects that can contribute to the intensities of each distinct MRI modality (13,14).
Studies have confirmed the reliability of the T 1 /T 2weighted MRI intensity ratio R as a trustworthy measure reflecting relative intracortical myelin content after TBI (15). Specifically, in TBI subjects, cortical maps of R are not only consistent with maps of myelin content inferred from histopathology, but also confirmatory of the fact that TBI is associated with relatively lower cortical myelin content in cross-section (15). To our knowledge, no detailed longitudinal study of post-traumatic cortical myelin change has been published.
Due to the potential to stimulate remyelination before axonal injury occurs, myelination-related therapies are being investigated to improve mTBI outcomes (16). Furthermore, better understanding of the effects of mTBI on brain architecture can improve surgical intervention (17,18). For these and other potential reasons, the study of posttraumatic demyelination is clinically relevant. This study quantifies how the T 1 /T 2 MRI intensity ratio R of the gray matter (GM) changes across the first ∼6 months after mTBI relative to age-and sex-matched healthy controls (HCs). Whereas other studies have mapped R after TBI, this is the first study that leverages MRI to quantify post-traumatic cortical myelin changes and to map the spatial patterns of these changes.

Participants and Data Acquisition
This study was conducted with the approval of the Institutional Review Board at the University of Southern California and in compliance with both the Declaration of Helsinki and the US Code of Federal Regulations (45 C.F.R. 46). Written informed consent was provided by all participants. A total of 80 HCs scanned at two timepoints were included from the Alzheimer's Disease Neuroimaging Initiative (ADNI; http://adni.loni.usc. edu). ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial MRI, positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of MCI and early AD. For up-todate information, see www.adni-info.org. Demographics of the HC cohort are listed in Table 1. The APOE genotypes were ε2/ε3 (N = 9 participants, i.e., 11% of the HC sample), ε2/ε4 (N = 1, Abbreviations: ADNI, Alzheimer's Disease Neuroimaging Initiative; APOE, apolipoprotein E; BTACT, brief test of adult cognition by telephone; df, degrees of freedom; DWI, diffusion weighted imaging; EVMD, episodic verbal memory delayed; EVMI, episodic verbal memory immediate; FS, FreeSurfer; GM, gray matter; ISI, interscan interval; IR, inductive reasoning; MPF, macromolecular protein fraction; MCI, mild cognitive impairment; mTBI, mild traumatic brain injury; PET, positron emission tomography; PS, processing speed; T E , echo time; T I , inversion time; T R , repetition time; VF, verbal fluency; WM, white matter; WMS, working memory span. i.e., 1%), ε3/ε3 (N = 47, i.e., 59%), ε3/ε4 (N = 22, i.e., 28%), and ε4/ε4 (N = 1, i.e., 1%). T 1 -weighted MRIs were acquired using a magnetization-prepared rapid acquisition gradient echo sequence (repetition time (T R ) = 2,400 ms; inversion time (T I ) = 1,000 ms; percentage sampling = 100; matrix size = 192 × 192; slice thickness = 1.2 mm). T 2 -weighted MRIs were acquired using a spin-echo sequence (T R = 3,000 ms; T E = 100 ms; slice thickness = 3.0 mm). A total of 70 HC subjects were scanned using a GE-manufactured scanner, while 10 were scanned using a Siemens-manufactured scanner. The myelin mapping procedure used here is valid only for examination of the cortex, whereas diffusion weighted imaging (DWI) is less accurate in GM than in white matter (WM) due to the relatively higher isotropy of water diffusion within GM (19). For this reason, DWI data were not included here. All participant data were deidentified and delinked. The mTBI group included 97 participants whose scans had been obtained both acutely (7.4 ± 3.2 days post-injury) and chronically (5.75 ± 0.3 months post-injury). Demographics of the mTBI cohort are listed in Table 1. For inclusion, mTBI volunteers had (a) one mTBI caused by a fall, (b) no clinical findings on acute T 1 /T 2 -weighted MRIs, (c) no clinical findings except cerebral microbleeds on susceptibility weighted imaging (a type of MRI on which hemorrhages and other iron-rich brain deposits are hypointense), (d) an acute Glasgow coma score greater than 12 upon initial medical evaluation, (e) loss of consciousness of fewer than 30 minutes, (f) post-traumatic amnesia of fewer than 24 hours, and (g) no clinical history of pretraumatic neurological disease or disorders, including dementia and MCI, psychiatric disorder and drug/alcohol abuse. T 1weighted MRIs were acquired using a magnetization-prepared rapid acquisition gradient echo sequence (repetition time (T R ) = 1,950 ms; echo time (T E ) = 3 ms; inversion time (T I ) = 900 ms; percentage sampling = 100; matrix size = 256 × 256; voxel size = 1 mm × 1 mm × 1 mm). T 2 -weighted images were acquired using a turbo spin-echo sequence (T R = 3,200 ms; T E = 222 ms; sampling percentage = 100; matrix size = 128 × 128; voxel size = 2 mm × 2 mm × 2 mm). All mTBI subjects were scanned using a Siemens-manufactured scanner. All participant data were deidentified and delinked. T 1 -and T 2 -weighted MRIs were acquired from each mTBI participant during two sessions held at least ∼5 months apart. mTBI subjects did not belong to the ADNI cohort.

MRI Preprocessing
Although calculating the ratio R of T 1 -to-T 2 image intensities facilitates the attenuation of receive field (B − 1 ) inhomogeneities that are inherent to T 1 and T 2 MRIs (12), bias field correction remains useful for alleviating remaining inhomogeneities and for improving the uniformity of MRI sensitivity to myelin content across image sequences (14,20). For these reasons, SPM12 software was used with default parameters to correct intensity inhomogeneity artifacts due to transmit fields, as detailed elsewhere (21). T 1 and T 2 volumes acquired during follow-up visits were registered rigidly to the T 1 and T 2 volumes acquired at baseline. The BRAINSFit algorithm (22) available in 3D Slicer (http://slicer.org) was used to apply successive coregistrations that improved the accuracy of the overall alignment procedure. These included (A) a rigid registration (6 degrees of freedom, df), (B) a rigid registration with scaling (7 df), (C) a rigid registration followed by scaling and shearing (10 df), (D) a fully affine registration (12 df), and (E) a B-spline registration with more than 27 df. To reduce registration errors, the transformation matrices associated with each registration were computed sequentially, multiplied across, and then applied to the MRI volumes as a single transformation performed in one step. By default, the percentage of image voxels sampled by BRAINSFit to determine the number and types of coregistrations is 0.2%. Here, this percentage was increased to 10% to improve precision.

MRI Segmentation
T 1 -and T 2 -weighted volumes were resliced to 1 mm 3 isotropic voxel resolution and segmented using FreeSurfer (FS) 6.0 software with default parameters, as detailed elsewhere (23,24). FS operations include (A) removal of non-brain tissue, (B) GM segmentation, (C) tessellation of the WM/GM boundary to generate triangular meshes of the meshes associated with the pial surface and with the GM/WM interface surface, and (D) surface topology correction. The quality of T 1 -weighted volume segmentations were enhanced by including T 2 -weighted volumes in the segmentation to improve pial surface reconstruction, the removal of tissue outside the brain, and the co-registration of volumes acquired using the two modalities. Regional cortical thickness and volumes were calculated by FreeSurfer and normalized using each subject's total intracranial volume. This normalization procedure removed the confound of head size, thereby enabling comparison across subjects.

Estimation of Relative Myelin Content
The ratio R of T 1 -to-T 2 -weighted volume intensities was calculated for each voxel after the volumes acquired at each wave had been co-registered within the same coordinate space. The FS GM mask was used to identify R values associated with GM voxels. The meshes representing the GM/WM interface surface and the pial surface were used for calculations that assigned an R value to each vertex on the GM/WM surface. The procedure resembled the one used by FS to calculate cortical thickness partly because correspondences between pial mesh vertices and WM mesh vertices were identified. These correspondences allowed us to compute, for each WM mesh vertex, the shortest distance between that vertex and the pial surface. Then, the line segment specifying the distance between the two surfaces was computed, subject to topological constraints to ensure proper pairing of vertices on the WM mesh to vertices on the pial surface mesh. For each such pair of vertices, the equation of the straight line that passed through the line segment connecting the two vertices was calculated. The value of the myelin measure R assigned to each vertex on the GM/WM surface was computed as a weighted average of the R values at the voxels traversed by this line segment. The averaging accounted for the number of GM voxels traversed and for the relative proportion of the line segment that lay within each of these voxels. For inter-subject averaging purposes, maps of R were co-registered to the FS average atlas.

Cognitive Measures Relative to R
Cognition was evaluated at approximately 2 weeks and 6 months post-injury in the mTBI cohort using the Brief Test of Adult Cognition by Telephone (BTACT) (25). The BTACT is a convenient phone-based cognitive assessment, comprising six classic measures of episodic verbal memory [both immediate (EVMI) and delayed (EVMD)], working memory span (WMS), processing speed (PS), verbal fluency (VF), and inductive reasoning (IR). Details on the administration, interpretation, and psychometric properties of these tests can be found elsewhere (26). Cognitive scores for the mTBI cohort are listed in Table 2.
Prior to analysis, mTBI cognitive scores x were converted to z-scores as z = (x − µ )/σ relative to the mean µ and standard deviation σ of x in a reference cohort comprising 4,513 healthy adults (2,027 males) aged 28-84 years (µ ± σ = 55.80 ± 12.31 years) who had participated in the Midlife in the United States (MIDUS) longitudinal study. Because BTACT cognitive scores were not normally distributed, Spearman's rank correlation was used to examine associations between R and cognition. We calculated correlations between (A) each measure of relative myelin content (i.e., R at each timepoint and the percentage change in R across timepoints) and (B) each cognitive variable assessed at every timepoint (EVMI, EVMD, WMS, VF, IR, and PS). Because BTACT scores are intercorrelated (27), a false discovery rate correction was implemented using the procedure of Benjamini & Yekutieli (28) for variables with dependencies.

Statistical Model
To compare R across groups (mTBI participants vs. HCs) and time points (acute baseline vs. chronic follow-up), a generalized linear mixed model with repeated measures was implemented using the fitlme function in MATLAB (MathWorks, Inc., Natick, MA). This model is of the form y = Xβ + Zu + ε, where y is a matrix containing the observations of the dependent (response) variable R, modeled as a continuous random variable with repeated measures. X is the design matrix of fixed effects and β contains the linear regression coefficients for the fixed effects. The three factors in the fixed model are sex (nominal variable), scanner manufacturer (discrete variable; 0 = Siemens; 1 = GE), and diagnosis (ordinal variable; 0 = absence of mTBI, i.e., HC; 1 = presence of mTBI). Age (continuous random variable with Gaussian distribution) and its higher powers up to fourth order are confounding variables. To identify the most appropriate fixed-effects model, the compare function in MATLAB was used to implement stepwise forward model selection. The compare function calculates the Akaike and Bayesian information criteria of two input models and uses a likelihood ratio test to determine which model best explains the underlying data without overfitting. Based on this strategy, the variables included in X were diagnosis, age, sex, scanner manufacturer, and the age × sex interaction. Because 10 HCs had been scanned on a Siemens-manufactured scanner, the column of X coding the scanner manufacturer was not colinear with that accounting for diagnosis effects.
In our model y = Xβ + Zu + ε, Z is the design matrix of random effect variables accounting for effects pertaining to interscan intervals (ISIs, i.e., by the time intervals t between baseline and follow-up scans). These random-effect variables are modeled by normally distributed, continuous random variables grouped by subject. In the HC group, the ISI models the amount of natural aging occurring between baseline and follow-up. In the mTBI group, it additionally models injury chronicity because baseline scans are acquired shortly after traumatic events. At baseline, t = 0 months in both groups; at follow-up, t varies across subjects. To account for ISI variability across subjects without regressing out the main effect of time, the deviation t i = u (t)− t i of the ISI t i of each subject i from the ISI grand mean u (t) was calculated first. These deviations were included in the generalized linear mixed model as random effects. The deviations from the ISI grand mean (rather than each subject's deviation from her/his group mean) were included because the two groups have different ISI means. For this reason, accounting for deviations from the individual group means (rather than from the grand mean) would have regressed out some of the interaction between the main effects of time and diagnosis. However, by accounting for deviations from the ISI grand mean, the confounding effect of ISI was treated identically across both groups within the generalized linear mixed model.
The vector u contains the random effects of t described above; the covariance matrix of u is G. The mathematical object ε is a vector of random errors with covariance V. Assuming normality, u ∼ N(0, G), ε ∼ N(0, V) and cov(u,ε) =0. Thus, our full model includes: (A) The within-subject response variable R with repeated measures that correlates with cortical myelin content across time. (B) The between-subject predictor variables coding for factors with fixed effects (i.e., sex, scanner manufacturer, and diagnosis). (C) One between-subject covariate (age at injury). (D) The between-subject interactions among fixed-effect predictor variables (e.g., sex × diagnosis, age × diagnosis, etc.). (E) The random-effect variables, grouped by subject, accounting for ISI.
The fixed-effect portion of our model includes an intercept for the grand mean, which is part of Xβ. The random-effect component Zu of the linear model organizes, by subject, the relationships between t, defined above, and R measurements. In this approach, the random-effect relationships are modeled by the linear equations y m = b 0m +b 1m t mp , where m is the index of each grouping variable (subject), and m = 1, ..., n since there are n subjects. The random-effect intercept b 0m and slope b 1m have prior distributions b 0m ∼ N (0, σ 2 0 ) and b 1m ∼ N (0, σ 2 1 ), respectively. The quantity t mp is the ISI deviation of subject m, at visit p, from the grand mean of t. The observation error term ε im has the distribution ε im ∼ N (0, σ 2 ), where i = 1, ..., np. The MATLAB fitlme function models the covariance matrix of a mixed linear system as a full matrix decomposed using Cholesky parametrization. It then uses maximum likelihood estimation to fit the unknown parameters using a trust region-based quasi-Newton optimizer with a relative tolerance of 1.0 × 10 −6 on the gradient of the objective function, an absolute tolerance of 1.0 × 10 −12 on step size, a maximum number of 10 5 iterations, and with an internally defined default value to start iterative optimization. To verify solution optimality, the algorithm checks the positive definiteness of the objective function's Hessian with respect to unconstrained parameters at convergence.

Hypothesis Testing
The probability density function of the average u(R), computed across both cortical vertices and subjects, was plotted for each diagnostic group (HC, mTBI) and timepoint. The average u(R) computed only across subjects was plotted on the cortical surface both before and after averaging across each cortical region. To quantify the size of within-subject effects, we first calculated R (the change in the relative myelin content R) as a percentage of R at baseline (R b ). Defining R f as the value of R at follow-up, it follows that R = (R f − R b )/R b . The calculation of R was implemented separately for HCs and for mTBI participants, while partialing out the effects of sex and age at injury. In other words, the adjusted values of R b and R f were obtained by regressing out the effects of both sex and age at injury. R was subsequently computed from the age-and sex-adjusted values of R b and R f . Of note, the random effects of ISIs were not partialed out when computing the adjusted values of R b and R f because doing so would have regressed out (much of) the random effect of time, which is of interest. Instead, to account for ISIs, the confoundadjusted value of R was converted to a percentage change in R observed across the mean ISI of the mTBI group. This procedure allowed us to quantify the change in R relative to a reference baseline group (HCs), independently of age at injury, sex, and of each subject's ISI. The probability density function of u( R) computed across both cortical vertices and subjects, was plotted for each diagnostic group. The average u( R) computed only across subjects was also plotted on the cortical surface.
To facilitate comparison of HCs and mTBI participants, the values of R and R were computed at each vertex of the cortical mesh, averaged across the participants in each diagnostic group, and then averaged across the cortical mesh vertices within each cortical parcel (29). The null hypothesis u(R TBI ) = u(R HC ) of no main effect of diagnosis on R was tested while partialing out the effects of all other confounding variables and their interactions. A paired-sample t test was used to test the null hypothesis u (R b ) = u(R f ) of no difference in R between baseline and follow-up within groups. This hypothesis was tested after regressing out the statistical effects of age at injury, sex, scanner manufacturer, and ISI. Welch's two-sided t test for the comparison of heteroskedastic samples was used to test the null hypothesis u( R TBI ) = u( R HC ) of no difference in R between HCs and mTBI participants independently of age at injury, sex, scanner manufacturer, and ISI. Because the distribution parameters of HCs' ISI differed from those of mTBI participants, the expected value assumed by R HC across the mean ISI of the mTBI participants was interpolated from the two repeated measures of R observed in HCs. This interpolation was based on a linear model which assumed a constant rate of change in R and which yielded the expected value of R HC across ∼6 months (i.e., the mean ISI of mTBI participants). Statistical hypotheses were tested for all cortical regions, averaged across hemispheres. All hypotheses were tested at the significance level α = 0.05. To control the family-wise error rate, corrections for multiple comparisons were implemented using a standard approach (30). The approach described above was also used for the analysis of regional changes in cortical thickness and volume.

Visualization and Comparison of Results
The difference in myelin percent loss between the two groups (mTBI and HC) is equal to u ( R TBI ) − u( R HC ), a quantity that was mapped on the cortex. To further explore these relative demyelination differences between groups, we computed the vulnerability v of each cortical region r relative to the set S of all other cortical regions as [u ( R r ) − u ( R S )] /u( R S ), where R r is R for cortical region r, and R S is R across all other cortical regions, i.e. across S. A positive value of v indicates that, on average, cortical region r experiences v% more demyelination than the rest of the cortex. Due to differences in scan parameters, direct comparison of R values between diagnostic groups may be problematic. Instead, to investigate the acute effects of mTBI on R within cortical parcel j, we computed the z-score z TBI = [R TBI − u (R TBI )]/σ (R TBI ) for TBI participants and similarly z HC for HCs. We then calculated the between-group difference in z-scores z = z TBI − z HC to compare mTBI participants' myelination to that of HCs. Using this procedure, we identified regions whose baseline R deviated excessively in mTBI subjects from its expected value in HCs. Advantageously, this approach does not assume that R has identical distribution parameters (u, σ ) within each group.

Age Sensitivity Analysis
Due to the difference in the age distributions of the groups, a sensitivity analysis was implemented to evaluate the regression model's effectiveness in partialing out age effects. All HCs younger than the oldest mTBI subject (i.e., younger than 79 y) and all mTBI subjects older than the youngest HC subject (i.e., older than 60 y) were selected to create subgroups whose participants whose ages had identical ranges. The HC subgroup consisted of 63 participants (age µ = 73 y, σ = 3 y) while the mTBI subgroup consisted of 17 participants (age µ = 71 y, σ = 4 y). The analysis described in subsection on hypothesis testing was repeated using only these subgroups.  Figure 1A displays, for all groups and timepoints, the probability density functions (PDFs), over cortical vertices, for the subject average of the confound-adjusted myelin ratio u(R). Figure 1B displays the PDFs of u( R), similarly computed. In mTBI participants, on average, u(R) decreases by 20.0% across the follow-up interval; this decrease is significantly greater than that of HCs across the entire cortex (all p < 0.05).   Pertaining to the analysis of acute mTBI effects described in subsection on visualization and comparison of results, on average, after accounting for age, sex, and other covariates, the largest z-score differences are localized to the occipital poles ( z = −0.24) and paracentral lobules and sulci ( z = −0.19, Figure 4). These regions, then, had lower u(R) values relative to u(R TBI ) than expected in the HCs. Figure 5 depicts cortical maps of u( R), according to the analysis described in the subsection on hypothesis testing. Whereas, both groups exhibit similar spatial patterns of cortical myelin loss, mTBI participants' demyelination is significantly more severe throughout the entire cortex (all p < 0.05). In HCs Figure 5A), regions with the smallest average decrease in u (R) include frontomarginal (1.9% loss on average), inferior occipital (1.86%), and transverse frontopolar (1.74%) gyri/sulci; middle occipital sulci (1.73%); and the occipital poles (0.98%). Regions with the largest average myelin loss include the temporal poles (3.88%); the middle-anterior (3.75%) and middle-posterior (3.60%) parts of the cingulate gyri and sulci; and the parahippocampal (3.45%), short insular (3.29%), and middle/superior temporal (3.28%) gyri. By contrast, in mTBI subjects (Figure 5B), myelin loss is more widespread and greater than in HCs. Across the followup period, all cortical regions are found to undergo decreases in u(R) larger than 10%. Structures with the greatest average myelin losses include the temporal poles (26.3%); the cingulate gyri/sulci (24.2%); and the orbital (23.9%), middle temporal (23.0%), angular (22.3%), inferior frontal (22.0%), superior frontal (21.9%), parahippocampal (21.4%), lingual (20.7%), and supramarginal (20.2%) gyri. These u( R) values are listed in Figure 6A. Whereas all cortical regions undergo significantly more demyelination in mTBI subjects compared to HCs, the extent of this phenomenon varies across regions. The smallest difference between groups involves myelin at the occipital poles (9.90% more demyelination, on average, in mTBI participants than in HCs). Figure 7, obtained from the analysis described in subsection on visualization and comparison of results, highlights regions that undergo the greatest myelin loss in mTBI subjects compared to HCs: the middle-anterior (22.7%) and ventroposterior (22.6%) parts of the cingulate gyri and sulci; the temporal poles (22.4%); the pericallosal sulci (21.6%); and the orbital (20.7%), short insular (20.3%), and superior temporal gyri (20.1%).

Regional Vulnerability to Demyelination and Atrophy
In mTBI participants, across the ∼6-month follow-up period, occipital regions experience 22.0% less demyelination, on average, than the rest of the cortex, whereas cingulate, temporal, and insular regions experience more demyelination, on average, than the rest of the cortex (19.8, 16.5, and 14.3%, respectively). In HCs, across the follow-up period, occipital regions experience 32.5% less demyelination than the rest of the cortex, whereas  Figure 1. Similarly, within each group, there is more myelin at the first timepoint compared to the second timepoint, even after adjusting for individual interscan intervals. However, whereas HCs experience relatively modest decreases in myelin with time, there is considerably more such loss in mTBI participants, particularly in dorsolateral frontal, lateral temporal, and especially occipital regions.
Frontiers in Neurology | www.frontiersin.org  Darker blue indicates that z < 0, i.e., that a cortical region has a lower z-score in mTBI participants than in HCs. For example, at baseline, the occipital poles and paracentral lobules and sulci exhibit substantially lower z-scores in mTBI subjects than in HCs. cingulate, temporal, and insular regions experience more demyelination than the rest of the cortex (27.7, 27.4, and 18.4%, respectively). In mTBI subjects, across the follow-up period, no change in mean regional cortical thickness or volume survived multiple comparison correction (q < 0.05). HCs experienced no significant differences in mean regional cortical thickness or volume, except for a significant decrease in cortical thickness within the inferior part of the left precentral sulcus. Figure 6B, illustrating the outcome of the analysis described in subsection on age sensitivity analysis, highlights that regional u( R) values computed across the full HC cohort are, overall, of similar magnitudes to those computed across the age-matched mTBI and HC subgroups. The ordering of regions is also similar, with the temporal and occipital poles consistently reported as being among the most and least demyelinated regions, respectively.

Correlation of R With Cognitive Measures
At acute baseline, the average of R across the cortex was negatively correlated with WMS (r = −0.223, p = 0.047) and with IR (r = −0.248, p = 0.027); however, these findings did not survive correction for multiple comparisons. No other significant correlations were observed at either timepoint, nor were any correlations between cognition and u( R TBI ). The chronic follow-up u(R) across the anterior transverse temporal gyrus was correlated with the WMS score at acute baseline (r = −0.311, p = 0.005). There was also a significant correlation between the chronic follow-up u(R) across the anterior temporal gyrus and the IR score recorded at acute baseline (r = −0.311, p = 0.005). A similarly significant correlation was found between the chronic follow-up u(R) across the anterior parts of the cingulate gyrus and sulcus and the IR score recorded at acute baseline (r = −0.308, p = 0.006). These last two findings did survive correction for multiple comparisons. The correlations of u (R) with cognition are listed in Table 3.

Interpretation and Implications
According to this study, independently of age and sex, mTBI participants exhibit significantly lower average R values at 6 months follow-up than at the acute timepoint, and significantly more demyelination than HCs. This suggests not only that mTBI can have lasting effects on myelin levels, but also that the size of such effects increases with time within the first 6 months after injury. Reassuringly, our HC myelin maps largely resemble those of previous studies (12,31). For example, unimodal sensory areas are highly myelinated whereas multimodal sensory and association areas are lightly myelinated. These trends are consistent with histological findings (31). In HCs, many regions with lower relative myelin content at baseline (e.g., temporal FIGURE 6 | The change R in the myelination ratio R for HCs (blue) and mTBI participants (orange), across both (A) the full cohort and (B) a sub-cohort containing only mTBI and HC subjects whose age ranges overlap. R was calculated as R = (Rf − R b) /R b , where R b is the R value at acute baseline and R f is the R value at chronic follow-up. R values were averaged across subjects and then across all cortical vertices within each region. Values are ordered in the ascending order of the full HC cohort's R values. Abbreviations: G, gyrus; HC, healthy control; mTBI, mild traumatic brain injury; S, sulcus. and frontal association areas) have been found to undergo relatively more demyelination than heavily myelinated regions (e.g., primary sensory areas), in agreement with findings from cross-sectional studies (20,32).
According to the z-score formula z = z TBI −z HC (subsection on visualization and comparison of results), a lower z-score in mTBI participants (i.e., z < 0) can reflect one of three scenarios: (A) R TBI < R HC and u (R TBI ) = u(R HC ); (B) R TBI = R HC and u (R TBI ) > u (R HC ); (C) both (A) and (B). The preponderance of the evidence suggests that myelin content decreases after mTBI (2,6). Therefore, when z < 0, u (R HC ) > u(R TBI ). At baseline, lateral occipital regions FIGURE 7 | Cortical map of group differences in R, corrected for multiple comparisons. These group differences convey deviations in mTBI-related demyelination from typical aging-related demyelination. Group differences were calculated as R TBI − R HC , where R TBI and R HC are R for the mTBI and HC groups, respectively. Darker blue shades indicate that mTBI participants experience a greater percentage of myelin loss than HCs. The former experience a greater relative amount of demyelination in temporal, cingulate, and insular regions, but also a relatively smaller amount in the pre-and postcentral gyri, as well as in occipital regions.
exhibit lower z-scores in mTBI participants than in HCs (i.e., z < 0, Figure 4) such that, at baseline, lateral occipital regions likely have lower R in mTBI subjects than in HCs. Conversely, in both HCs and mTBI participants, these occipital areas were also found to be among the least vulnerable to demyelination based on the percentage difference in myelin between timepoints in each region (Figure 5). At baseline, temporal, cingulate, and insular regions have higher z-scores in mTBI subjects than in HCs, indicating that these areas may be among the least vulnerable to acute demyelination. However, these regions are also among the most vulnerable to chronic demyelination. Thus, our results indicate that cortical regions can differ significantly in their vulnerability to acute vs. chronic mTBI. The differential vulnerability of brain regions to demyelination in acute vs. chronic mTBI suggests that the mechanisms underlying acute post-traumatic demyelination differ, at least in part, from those resulting in chronic demyelination.
Although demyelination is significantly more severe in chronic mTBI than in typical aging, the spatial pattern of regions most or least vulnerable to demyelination overlaps substantially across these groups. Such similarities may indicate that, typically, these brain regions are similarly vulnerable to both aging-and mTBI-related demyelination mechanisms, or that demyelination is caused by similar factors in both groups. This supports the hypothesis that chronic mTBI-related demyelination is linked to secondary mechanisms of brain damage, which typically include neuroinflammation and excitotoxicity that are also manifest in typical aging (33)(34)(35). Furthermore, the regions most vulnerable to demyelination across both groups (e.g., insular, temporal, and cingulate cortices) overlap greatly with those that undergo significantly more demyelination after mTBI than is typical (Figure 7). Conversely, the regions least vulnerable to demyelination across both groups (e.g., occipital, frontomarginal) undergo significantly less demyelination than the average amount expected after mTBI. This suggests the testable hypothesis that chronic mTBI effects on myelin may compound non-linearly in these areas, possibly due to interplay between various secondary injury mechanisms.
The paucity of significant correlations between cognitive scores and R values is likely influenced by cognitive recovery (36,37). Post-traumatic recovery of myelin occurs to some degree in WM (9) but is understudied in the GM. Primary injury mechanisms can result in axonal damage and cell death before remyelination can begin (38), and studies of intracortical remyelination in multiple sclerosis reveal that, even when myelin recovery occurs, the patterns of myelination may be altered pathologically (39). The negative associations between IR, WMS, and R are surprising, but likely reflect that IR and WMS are not exclusively dependent on focal myelin profiles. IR and WMS are cortically distributed processes that require communication across multiple GM areas via WM connections. Thus, it is possible that R reflects short-range local communication better than long-distance cortico-cortical processing of stimuli or neurally encoded information.
Occasionally, the spatial pattern of mTBI-related demyelination mapped here deviates from that observed in typical aging-related demyelination, e.g., in the case of sensorimotor cortex. Pre-and postcentral gyri (which harbor the primary motor and somatosensory cortices, respectively) experience substantially more typical aging-related demyelination, but substantially less mTBI-related demyelination than most other regions. In this respect, there are notable similarities between mTBI-and AD-related demyelination patterns. Although the latter remains understudied in the cortex, it is known that white matter tracts linking temporal areas to the rest of the cortex undergo significantly more demyelination in AD, whereas tracts proximal to motor and somatosensory cortex are relatively spared (40). Furthermore, AD causes excessive cortical thinning of orbitofrontal, anterior temporal, parietal, cingulate, and insular regions (41)(42)(43), while sparing somatosensory regions (44). This cortical thinning pattern mirrors the demyelination pattern documented here in chronic mTBI ( Figure 5B). AD-related cortical thinning is partly due to neuroinflammation of the cortex (45), which is known to be both a cause and a consequence of demyelination, across both AD and TBI (33,46). In this context, our findings suggest that there is substantial overlap between the brain's spatial patterns of vulnerability to (A) AD-related cortical thinning, (B) AD-related white matter demyelination, and (C) mTBI-related intracortical demyelination. Thus, because mTBI increases the risk of AD and related dementias (47, 48) through uncertain mechanisms, further study on how post-traumatic demyelination patterns can translate into higher AD risk is warranted. It is possible that edema lowers R after mTBI. Edema may cause hyperintensity on T 2 -weighted MRI, leading to lower R values, and both vasogenic as well as cytotoxic edema can become manifest after mTBI (49). Cortical thickening can be a correlate of post-traumatic edema, because this latter phenomenon is often due to acute inflammation leading to intra-and extracellular fluid accumulation that translates into larger regional volumes and thicker cortex despite frequent loss of neurons (50). Given that the mTBI subjects in this study experienced no significant change in cortical thickness or volume across the follow-up period, edema may not have played a substantial role in causing reductions in R across ISIs.
It has been argued that R may correlate with neuroanatomical features other than myelin. Correlations between R and dendrite density, axon caliber, and cell protein size have been reported, as has a negative relationship between R and mitochondrion gene variants (51,52). The latter association may imply abnormally high mitochondrial energetics in brain regions with low R in this study (e.g., parts of the temporal, cingulate, and insular cortex). Abnormally high mitochondrial energetics are known to contribute to neurotoxicity (53). On the other hand, the association between protein size and R, which is partially related to the negative correlation between R and proteasome-expressed genes (52), could indicate malignant proteolysis, which has been described as a primary effect of TBI that accompanies traumatic axonal injury (54). Thus, R may correlate with markers of both primary and secondary injury after mTBI. Future studies should explore complementary measurements that can refine the interpretation of R values and their changes across the brain.
Whether based on the full cohort or on the age-overlap subgroups, our within-group R comparison is qualitatively and quantitatively consistent. This suggests that our regression model was effective in partialing out the statistical effect of age, despite our use of a relatively younger mTBI cohort and of a relatively older HC cohort. This sensitivity analysis also suggests that age distribution discrepancies between HCs and mTBI participants did not have a substantial effect on our ulterior analysis.

Comparison With Other Studies
The relative amount of post-traumatic demyelination reported in this study is noteworthy. Although diffusion MRI studies of mTBI do not typically report brain circuitry changes of similar magnitudes (55,56), such studies focus on WM changes rather than on GM changes due to the far more isotropic diffusion of water within GM compared to WM (19). Furthermore, the extent to which diffusion MRI measures of WM integrity are related to myelin content in the GM is poorly understood, thus precluding interpretation of our findings relative to diffusion MRI studies. Magnetization transfer ratio studies similarly tend to find smaller effects than ours, but these studies are also focused primarily on WM (57). One of the few imaging other measures to survey intracortical myelin is the macromolecular protein fraction (MPF). To our knowledge, only one study has used MPF to investigate myelin content after human mTBI. That study reported MPF reductions of 18-36% in precentral, anterior cingulate, medial orbital, subcallosal, lingual, superior frontal, and inferior frontal gyri (58). However, these reductions were observed across periods much longer than the typical ISI in our study (i.e., ∼10 years versus ∼6 months), and in persons who had sustained several mTBIs rather than just one, as in our case. Thus, interpretation of our results relative to available MPF findings requires further studies.
It is illustrative to compare the findings in our HC cohort to those in a cross-sectional lifespan study of R by Grydeland et al., who studied an HC cohort different from ours (20). Despite similar findings in terms of regional patterns of demyelination (e.g., greater demyelination in cingulate/temporal cortex, less demyelination in sensory cortex), Grydeland et al. reported more gradual declines in R than are seen in our study at certain cortical regions. Discrepancies between the two studies may be due to methodological differences, such as our inclusion of bias field correction and our mapping of demyelination in a longitudinal cohort rather than crosssectionally.
Importantly, we must remind the reader that the parameters of the mathematical function describing the relationship between R and true cortical myelin content are not known with high precision. For this reason, we cannot determine whether our computed values of R reflect the true percentage change in cortical myelin that would be determined by direct measurement. It is possible, for example, that R scales linearly with true myelin content, or that this scaling is non-linear. However, because the precise coefficients of either model are unknown to us, our reported R values may convey a measure of relative demyelination extent rather than true values of myelin content change.
Only one other study known to us (15) calculated R in persons with mTBI, finding a significant negative association between R and the number of mTBIs sustained by participants. The negative association was strongest in lateral occipital areas, which is inconsistent with our findings on regional vulnerability to chronic mTBI (Figure 5B), but consistent with those on acute mTBI (Figure 4). This discrepancy may imply a shared mechanism of demyelination between repeat and acute mTBI, but further research is needed to clarify the matter. One interpretation is that primary injury mechanisms drive acute post-traumatic demyelination and are upregulated in repeat mTBI because of multiple distinct injuries (55). Future studies should integrate T 1 -and T 2 -weighted MRI with other approaches that can distinguish primary from secondary injury effects upon the cortex.

Limitations
HCs and mTBI participants were not imaged as part of the same study. Thus, on average, HCs are typically older than mTBI participants and their ISIs are longer as well. Although these group differences were accounted for in the statistical analysis, we acknowledge that it would have been preferable for the HC and mTBI groups to be better matched according to these variables. In our HCs, 30% had at least one ε4 allele, compared to ∼13% in the general population (59). The ε4 genotype is known to be associated with higher AD risk and more severe demyelination (60). Thus, differences in genotype between cohorts could partly explain some group differences reported here, particularly if the mTBI group had fewer ε4 allele carriers than the HC group. Because the APOE genotype of mTBI participants was unknown to us, its statistical effect could not be accounted for. We acknowledge this as a limitation of our analysis; future studies should quantify whether and how APOE genotype affects posttraumatic demyelination.
The drawbacks of estimating myelin content from the ratio R of T 1 /T 2 -weighted MRIs should be acknowledged. Although R is significantly correlated with myelin-related gene expression in the cortex, R is correlated more strongly with axon caliber and oligodendrocyte markers (52) and the myelin water fraction index has been proposed by Arshad et al. (61) as a preferable index of myelin content, although these authors compared R and myelin water fraction only from the standpoint of their abilities to quantify myelin in the white (rather than gray) matter. A later study (62) found that myelin water fraction correlates poorly with R in white matter, but more strongly in gray matter, such that our use of R in this study of cortical myelin is likely justified. Aside from these considerations, however, we acknowledge that estimating myelin content from the ratio of T 1 /T 2 -weighted MRIs can be affected by intensity scaling disparities across individuals, scanners (61), and possibly even across neurological conditions (63), although evidence for this latter effect is equivocal (64). Specifically, Pelkmans et al. found higher R values in AD compared to cognitively normal persons. However, as they discuss, these results were likely influenced by the iron deposition that accompanies AD. Iron is reflected as hyperintensity in T 1 images, resulting in an increased R value. mTBI, on the other hand, likely does not cause the same level of cortical iron deposition within 6 months after injury. Acute iron deposition in mTBI is mostly heme-bound, manifesting as cerebral microbleeds (65), which are typically small and few (66), although their incidence increases with age (67). Non-heme iron accumulation after mTBI is understudied in humans but has, thus far, been found only near ventricles, not in cortical GM (65,68).
Because follow-up scans were registered to baseline T 1weighted scans for each subject, the follow-up T 1 -weighted scans were subjected to more interpolation than baseline T 1weighted scans. Due to the high sampling percentage used by the registration procedure and to the high fidelity of the registration between baseline and follow-up scans, this discrepancy is likely to have had very limited effect upon MRI intensities. However, we acknowledge that applying an identical interpolation to baseline T 1 -weighted scans may have been preferable. Although mTBI and HC participants were not scanned as part of the same study, the scanner parameters were similar across mTBI and HC groups. Our previous studies (66, 69) on a sample that overlapped substantially with that in the current one quantified differences in aging between mTBI and HC groups and found no significant differences between these groups pertaining to the statistical effects of sequence parameters and scanner types. Additionally, comparisons of R were unlikely to be impacted by this confound, as scanner parameters were identical across follow-up for all subjects. Nevertheless, we acknowledge that scanning both mTBI and HC participants using the same scanner and sequence parameters would have been preferable. The z-score analysis implemented for the baseline comparison circumvented the issue by comparing the distribution of R values between groups, but we acknowledge that this analysis relied on the assumption that the effects of scanner parameters on R were mainly consistent across the cortex. If scanner parameters had significant variation in their effects on R at different cortical regions, that could potentially complicate the interpretation of acute mTBI effects on relative myelin content. This confound is expected to have had minimal impact, especially due to the nature of the T 1 /T 2weighted ratio, which increases contrast for myelin specifically, thereby reducing the variation in response to other tissue features that T 1 -or T 2 -weighted imaging may exhibit alone (12,14). The regression model partialed out the statistical effect of withingroup variation related to scanner differences, a confound that only impacted HCs. The cortical patterns of R was very similar across HCs scanned using GE and Siemens scanners, and there was no significant difference (p < 0.05) in either µ(R) or σ (R) across scanner manufacturers. The sequence parameters were consistent across scanner manufacturers.

CONCLUSION
Mapping the spatiotemporal profile of demyelination after mTBI is clinically relevant. This study suggests that chronic mTBIrelated intracortical demyelination has a distinct spatial profile from that observed in the acute phase of mTBI. Furthermore, post-traumatic demyelination is substantially more severe than typical aging-related demyelination. These and related findings of our study have implications for elucidating the mechanisms whereby mTBI-related demyelination occurs. Because posttraumatic cortical demyelination can affect clinical outcomes and may even reflect the risk for neurodegenerative diseases like AD, future studies should further investigate this phenomenon relative to other factors like age, sex, APOE genotype, cognitive function, and other biological factors. Future research should also investigate the association between our MRI measures and independent biomarkers of demyelination, e.g., serum levels of neurofilament light chain protein.

DATA AVAILABILITY STATEMENT
Primary data generated and/or analyzed during the current study are available subject to a data transfer agreement. At the request of some participants, their written permission is additionally required in some cases.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board at the University of Southern California. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
SM: study concept and design, algorithm implementation, data acquisition and analysis, and contributions to manuscript. NC: study design, data analysis and interpretation, and contributions to manuscript. VN: algorithm implementation and data acquisition and analysis. PI: data analysis and contributions to manuscript. AI: study design, manuscript revision and review, data analysis and interpretation, and implementation of statistical analysis. All authors contributed to the article and approved the submitted version. FUNDING grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California. The funders were not involved in the study design, collection, analysis, interpretation of data, the writing of this article or in the decision to submit it for publication.