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-to-date 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, 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₁-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₂-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₁/T₂-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 pre-traumatic neurological disease or disorders, including dementia and MCI, psychiatric disorder and drug/alcohol abuse. T₁-weighted 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₂-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₁- and T₂-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.

Although calculating the ratio R of T₁-to-T₂ image intensities facilitates the attenuation of receive field (B1-) inhomogeneities that are inherent to T₁ and T₂ 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₁ and T₂ volumes acquired during follow-up visits were registered rigidly to the T₁ and T₂ volumes acquired at baseline. The BRAINSFit algorithm (22) available in 3D Slicer (http://slicer.org) was used to apply successive co-registrations 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 co-registrations is 0.2%. Here, this percentage was increased to 10% to improve precision.

T₁- and T₂-weighted volumes were resliced to 1 mm³ 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₁-weighted volume segmentations were enhanced by including T₂-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.

The ratio R of T₁-to-T₂-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.

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.

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:

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 b0m ~ N (0, σ02) and b1m ~ N (0, σ12), 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⁻⁶ on the gradient of the objective function, an absolute tolerance of 1.0 × 10⁻¹² on step size, a maximum number of 10⁵ 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.

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 confound-adjusted 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.

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.

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.

Figures 1–3 illustrate the results obtained from the analysis described in the first paragraph of the subsection on hypothesis testing. 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).

Figure 2 displays cortical maps of R in HCs (2A: baseline; 2C: follow-up) and mTBI participants (2B: baseline; 2D: follow-up). Overall, cortical R patterns are similar across both groups and time points. The observed patterns of R bear similarity to those of Glasser et al. (12), who used a similar quantitative approach to map myelin content. For example, in HCs at baseline, u(R) is highest within the inferior occipital (u(R) = 1.25) and anterior transverse temporal (u(R) = 1.23) gyri, at the occipital poles (u(R) = 1.21), within precentral (u(R) = 1.16) and postcentral gyri (u(R) = 1.15), and within frontomarginal (u(R) = 1.15), subcentral (u(R) = 1.18), and transverse frontopolar gyri and sulci (u(R) = 1.19). Figure 3 displays HCs' cortical map of u(R) at baseline, using a color scheme resembling that of Glasser et al. In both groups and across both timepoints, u(R) values are relatively lower within most cingulate areas, except for the moderately myelinated dorso-posterior part of the cingulate gyrus (u(R) = 1.08). u(R) values are typically lowest in the ventro-posterior (u(R) and middle-anterior (u(R) = 1.01) parts of cingulate cortex. As expected, the insula has relatively low u(R), as do frontal and temporal lobes. The temporal pole is very lightly myelinated (u(R) = 0.96).

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 follow-up 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%).

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 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.

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.

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 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 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–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–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₂-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.

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 cross-sectionally.

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₁- and T₂-weighted MRI with other approaches that can distinguish primary from secondary injury effects upon the cortex.

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 post-traumatic demyelination.

The drawbacks of estimating myelin content from the ratio R of T₁/T₂-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₁/T₂-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₁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₁-weighted scans for each subject, the follow-up T₁-weighted scans were subjected to more interpolation than baseline T₁-weighted 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₁-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₁/T₂-weighted ratio, which increases contrast for myelin specifically, thereby reducing the variation in response to other tissue features that T₁- or T₂-weighted imaging may exhibit alone (12, 14). The regression model partialed out the statistical effect of within-group 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.