Aberrant levels of cortical myelin distinguish individuals with unipolar depression from healthy controls

1 INTRODUCTION

Unipolar depression (UD) is a leading cause of disability worldwide¹, with an economic burden of $210 billion dollars in the United States alone². Despite its impact, treatments for the disorder remain ineffective for many patients³. Thus, there is a pressing need to understand the neurobiological etiology of UD to facilitate the development of improved treatments and prevention strategies.

Unipolar depression is characterized by dysfunctional affective and cognitive processing⁴, including reduced executive functioning⁵, biased emotional processing⁶, and impaired reward processing⁷. Correspondingly, individuals with UD show aberrant activation during tasks which recruit these processes, including activation in the striatum, hippocampus, amygdala, orbitofrontal cortex, prefrontal cortex, insula, cingulate, and occipital cortex^8–11. In addition, a growing body of literature has reported structural abnormalities associated with depression in many of these same regions described above, both in grey matter^12–14 and in white matter^15,16. Meta-analyses of diffusion weighted imaging (DWI) studies have repeatedly found evidence for lower fractional anisotropy (FA) in depressed populations^17–19. More recently, studies using very large samples (i.e., the UK Biobank), as well as meta-analyses combining both published and unpublished data (i.e., the ENIGMA consortium), have observed widespread and replicable reductions in FA^15,16. Notably, white matter integrity in identified regions has been shown to correlate with the cognitive processes disrupted in depression, including processing speed^20,21, emotion regulation^22,23, and reward learning²⁴.

Emerging evidence suggests that cortical myelin may be impacted in individuals with UD and that it may partially mediate some of the cognitive processes that are impaired in affected individuals. For example, studies of post mortem brain tissue from donors with UD have observed reduced myelination, a reduction in the number of oligodendrocytes and other glia (cells whose functions include generating and maintaining myelin²⁵), and reduced expression of oligodendrocyte lineage genes^26,27. A study of individuals with treatment resistant depression revealed a reduced magnetization transfer ratio (MTR), which is thought to reflect lower myelin levels, in the cingulate cortex and insula²⁸. A recent study using R1 (1/T1) as a measure of myelination observed reduced whole-brain myelin, but no significant difference in cortical myelin in a handful of a priori bilateral regions between individuals with depression and healthy controls (HC)²⁹,

Developments in magnetic resonance imaging (MRI) methodology permit the examination of cortical myelin via the T1w/T2w ratio^30,31. Studies in population-based samples using this metric have found that lower myelin in the cingulate, orbitofrontal cortex, and middle temporal cortex correlated with poor sleep quality³², lower frontal-pole myelin and greater myelin in the occipital cortex correlated with neuroticism³³, and lower myelin in the motor and higher myelin in the insular, cingulate, prefrontal, and superior parietal cortices correlated with trait anxiety³⁴. While poor sleep, high neuroticism, and trait anxiety might represent concurrent symptoms of depression, prior studies have not systematically examined cortical myelin in participants with depression as compared to HC.

The goals of the present study were (1) to ascertain whether cortical myelin content (characterized by the T1w/T2w ratio) is predictive of unipolar depression (UD), and (2) to characterize the brain regions that are predictive of case/control status. Based on the prior studies mentioned above, we hypothesized that cortical myelin levels would distinguish individuals with UD from HC and that these differences will be especially pronounced in the prefrontal cortical (PFC), cingulate, parietal and occipital regions that support reward and emotional processing, which are dysregulated in UD^11,35.

2 METHODS

2.4 Data analyses

2.4.1 Clinical data analysis

HC and UD groups were compared on demographic and clinical variables using t-tests and chi-square tests. All analyses were conducted in R (https://www.r-project.org/).

2.4.2 Neuroimaging data processing

2.4.2.1 Subject-level preprocessing

Data quality was examined using mriqc version 0.15.1⁴⁶ and visually inspected. Participants were excluded from the analyses if visual inspection identified the presence of gross motion or scanner-related artifacts in the T1w or T2w images, or in the mriqc background noise images, (e.g., ghosting, blurring, ringing, banding, etc.). The distribution of mriqc image quality metrics (IQMs) reflecting scan noise in our study were compared to the distribution of IQMs collected by the mriqc server from other studies and available via its API, from 1046 T1w and 619 T2w images collected using similar parameters⁴⁷. Scans with IQMs beyond the interquartile range of the mriqc API data (median +/-1.5 x 75% quartile – 25% quartile) were flagged as potential outliers and were re-inspected. The distributions of the mriqc IQMs for the remaining participants were not different from the distribution of IQMs from the mriqc API (Supplemental Figures 1&2).

Each participant’s cortical myelin was characterized with the T1w/T2w ratio^30,48,49 using the PreFreeSurfer, FreeSurfer, and PostFreeSurfer minimal preprocessing pipelines for the human connectome project⁴⁸. Workbench v1.4.2 and HCPpipelines-4.1.3 were installed system-wide on a workstation with GNU/Linux Debian 10 operating system. The spin echo field maps collected in AP and PA phase encoding directions were used for bias field correction in PreFreeSurfer. Registration to standard space was achieved via MSMSulc⁵⁰ in PostFreeSurfer. FreeSurfer images and T1w/T2w myelin maps were visually inspected for errors and artifacts and images with gross errors (e.g., large regions of apparent low myelin in the occipital cortex due to the transverse sinus interfering with accurate identification of the pial surface) were removed from analyses. The resulting myelin maps were parcellated in Workbench using the 360 region Glasser Atlas⁴⁹.

2.4.2.2 Group-level preprocessing

Data quality in some brain regions is worse than in the other regions due to susceptibility artifacts⁵¹. This inconsistent quality of data may increase variability in myelin values across participants. The distribution of myelin values across all participants in each of the 360 regions was examined to identify regions whose variability was an outlier relative to variability in other regions. The coefficient of variation (sd/|mean|) was used to summarize the variability within each region. Rosner’s test for outliers^52,53 identified 11 outlier parcels with excessively high variation (Supplemental Table 1) including parcels located in bilateral hippocampus, entorhinal cortex, presubiculum, piriform cortex, and posterior orbitofrontal cortex complex, and the right subgenual cingulate (bilateral H, EC, PreS, Pir, pOFC, and right 25).These regions, which are known to suffer from an excess of susceptibility artifacts,^54,55 were removed from the analyses, leaving 349 parcels.

2.4.3 Neuroimaging data analysis

2.4.3.1 Elastic net and linear discriminant analysis

Neuroimaging data decoding studies can capitalize on complex relationships between variables, but their large numbers can present a challenge in deriving a generalizable model. The elastic net approach has emerged as a flexible tool for use with neuroimaging data⁵⁶ as it is able to reduce the influence of overly large coefficients and reducing the number of variables while generating multivariate models predictive of complex behaviors^57–59. Elastic net is a regularized regression which combines lasso and ridge regression (i.e., L1- and L2- norm regularization)⁶⁰. Ridge regression penalizes overly large coefficients, while lasso regression removes variables with small coefficients. Elastic net has two parameters: alpha (α) controls the balance between the ridge and lasso regularizations, and lambda (λ) controls the strength of regularization.

We used logistic elastic net regularized regression implemented in the R package glmnet⁶¹ to select variables that were most predictive of case/control status. To provide equal contribution of each penalty to the loss function, we used α=0.5. To avoid model overfitting and bias, we implemented nested cross-validation to identify the optimal λ parameter. A linear discriminant analysis (LDA)⁶² model was subsequently trained using selected variables to make out-of-sample prediction on held-out participants. This strategy is illustrated in Figure 1A and described in detail in Supplemental Methods. For each repetition of the nested cross-validation loop, two participants (1 UD and 1 HC) were held out. The rest of the sample was used to identify the optimal λ parameter which were then used to fit the elastic net model and select variables whose myelin levels were predictive of UD/HC status. These variables were then used to train an LDA model, which was tested on held-out participants. Results (model fit, variable selection, prediction accuracy) were evaluated for statistical significance using a permutation analysis (Figure 1B and Supplemental Methods).

• Download figure
• Open in new tab

Figure 1. Diagram of analysis steps

Conceptual depiction of analysis steps including: (1) a unique pair of one UD and one HC participant is held-out; (2) an elastic net regression is used to select variables; (3) the retained variables are used an LDA model predicting case/control status; (4) the LDA model is tested on the held-out sample; (5) this process is repeated for each of the n=1833 pairs of subjects; (6) for each held-out pair, the training procedure is repeated with 100 unique permutations.

2.4.3.2 Post-hoc analyses

2.4.3.2.1 The relationship between true-label and permuted-label sample demographics

To ensure that the permuted results were not due to changes in the internal structure of the permuted samples in terms of demographic variables (i.e., age, sex and IQ), we compared the age, sex and IQ values in the permuted samples with that in the true-label samples.

2.4.3.2.2 Association of LDA accuracy and cortical myelin with demographic and clinical characteristics

To further characterize the parcels selected by the logistic elastic net regression, we tested the association of participant group (i.e., UD vs HC) with cortical myelin. These analyses controlled for age, sex, and IQ. To assess the potential influence of confounding variables on model accuracy, we tested whether demographic (i.e., age, sex, and IQ) or clinical characteristics (i.e., HRSD-25 and MOODS-SR) were predictive of classification accuracy. This was done via a two-way ANOVA, which tested whether group (i.e., UD vs. HC), the variable in question, or their interaction, was associated with participant-wise accuracy⁶³. Additionally, a 1-way ANOVA tested whether clinical characteristics of UD participants (e.g., medication type, medication load, duration of illness) were associated with model accuracy only in UD participants. To assess the potential influence of confounding variables on variable selection, ANCOVAs tested the association of demographic and clinical variables with cortical myelin. These analyses controlled for age, sex, and IQ, and were performed across both UD and HC participants (i.e., 2-way Group x Variable ANCOVA), as well as only in UD participants (i.e., 1-way ANCOVA). Results for each variable were separately corrected for multiple comparisons using false discovery rate (FDR).

2.4.3.3 Exploratory analysis of a HC participant who was diagnosed with major depressive disorder 12 months from the baseline scan

One participant entered the study as a HC but was diagnosed with major depressive and generalized anxiety disorders sometime between 6 and 12 months after study onset. At the study visit at 12 months the participant had mild depressive symptoms. While this participant was excluded from all primary analyses described above, we thought it would be informative to conduct exploratory analyses investigating whether the myelin was predictive of the participant’s conversion from HC to UD. This exploratory analysis used the primary 86 participants and the variables selected in primary analyses (cortical myelin in 33 parcels and IQ) to train an LDA model. UD/HC status was then predicted at baseline (12 months prior to conversion) and at the 6-month follow-up (6 months prior to conversion).

3 Results

3.2 UD vs. HC LDA nested cross-validation classification accuracy

Cortical myelin levels and IQ distinguished UD from HC in subjects held-out during cross-validation, with an average accuracy of 68% (Figure 2A; sensitivity (UD): 63.8%, specificity (HC): 71.5%). The mean participant-wise accuracy across 84 classification loops ranged from 0% to 100%. Greater than 50% accuracy was achieved in 63 (73.6%) participants (26 (66.7%) UD, 37 (78.7%) HC, exact binomial test p<0.001). When demographic variables were excluded from the LDA classification, nested cross-validation achieved 69% accuracy (UD: 65.2%, HC: 73%). Notably, while the mean classification accuracy was higher in the HC group, this difference was not statistically significant (t=0.88, p=0.38), and similarly the proportion of participants classified at greater than 50% accuracy did not differ statistically between the two groups (χ²=1.03, p=0.31).

• Download figure
• Open in new tab

Figure 2. Classification accuracy and variable selection

A) Average subject-wise nested cross-validation classification accuracy in both healthy controls (HC) and participants with unipolar depression (UD). B) Variable selection frequency (percent of models where a given parcel was retained) in the true data (light blue) and permutations (green). The n=90 parcels which were retained at least once are shown. IQ, the only demographic variable selected, is labeled. Solid black lines represent the median of the variable selection frequency for all n=350 variables +/-3.5 x the interquartile range (IQR). See Supplemental Table 2 for the full list of variable selection frequencies.

In permutation analyses, no variable was selected in 66.9% of 183300 models. Within each participant, the proportion of models that did not select any variables ranged from 65.4% to 68.7%. Excluding instances when no variables were selected, the average participant-wise LDA accuracy in permutation analyses ranged from 41.2% to 59.2%, with an average accuracy of 50.5% (i.e., chance level). The test of whether age, sex and IQ in permuted samples differed from that in the sample with true labels showed that the demographics of the permuted groups differed (p<0.05, uncorrected) from the true-label groups only in 1.7% of cases.

3.3 Elastic net variable selection

True-label nested cross-validated elastic net models predicting diagnostic status (UD vs. HC) with 349 myelin parcels, age, sex and IQ, selected 90 myelin parcels (Supplemental Table 1) and IQ in at least one model (Figure 2B). Models selected between 9 and 68 variables, with a median of 17 variables. Permuted-label nested cross-validated elastic net models selected all predictor variables at least once, but no variable was selected by more than 3% of the models. In addition, 66.9% of models with permuted labels did not select any variable at all. Within permuted-label models where at least one variable was selected, the number of selected variables ranged from 1 to 113, with a median of 17 variables selected.

Given that the variable selection frequency in the models with permuted labels likely represents noise, we applied the criterion of the median + 3.5*IQR across all permuted variables (3.77%) as the cutoff value to separate the potential noise variables from ‘signal’ variables in the nested cross-validated elastic net with true labels. This latter analysis identified 33 out of 90 myelin parcels, plus IQ, across the true-label nested cross-validated elastic net models that were above the cutoff line (Table 2, Figure 2B). These 33 parcels included multiple regions in the orbitofrontal cortex, insula, cingulate, and frontal operculum, as well as regions in the auditory and visual cortices (Supplemental Figure 3). Of these parcels, six parcels (left 11l, left 7PC, left MST, left p10p, right FOP2, and right FFC) were retained in all 1833 models (100%), an additional 7 parcels were retained in more than 90% of models, 13 more parcels were retained in more than 10% of models, and 7 parcels were retained in less than 10% of models. IQ was retained by 85.8% of models. No model retained the age and sex variables.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2.

Association of unipolar depression with cortical myelin in selected parcels

Post-hoc analyses tested the association of diagnostic status with cortical myelin levels in the 33 selected myelin parcels (Table 2, Figure 3). These parcels include both regions where UD participants have lower mean cortical myelin than HC and regions where UD participants have greater mean cortical myelin, as well as regions where the two groups do not differ in their average level of cortical myelin. After FDR-correction for multiple comparisons, fourteen of the parcels showed evidence of significant differences between UD and HC participants (Table 2). These parcels included left 7PC, right TE2p, and left MST, where participants with UD had greater myelin levels, and left s32, left FOP4, right FFC, left RI, right 7AL, right FOP2, left 24dv, left p10p, left LBelt, left a24pr, and left 11l where participants with UD had lower myelin levels. (see Table 2 for an expanded description of each parcel). An additional 8 parcels showed nominally significant differences between UD and HC participants (Table 2; p<0.05 uncorrected), including left STSdp, left d23ab, left 6a, and right LBelt, where participants with UD had lower myelin levels, and left p32, right PH, left MT, and right RSC, where participants with UD had greater myelin levels. Parcels that showed a greater absolute mean difference between the groups were selected more frequently in the nested cross-validation analysis (r=0.7, p=4×10⁻⁶).

• Download figure
• Open in new tab

Figure 3. Association of Unipolar Depression with cortical myelin in selected parcels

A) The t-statistic for the association of unipolar depression (UD) with cortical myelin in the n=33 selected parcels. Regions where UD participants had greater average myelin than HC participants are in blue, and regions where HC participants had greater average myelin than UD participants are in red. B) The percent of models in which each parcel was retained. Regions are individually labeled (note that FFC is not visible, as it lies on the ventral surface).

4 DISCUSSION

Cortical myelin distinguished healthy control participants (HC) from those who had been diagnosed with unipolar depression (UD) with 68% accuracy in 2 subject hold-out cross-validations. Classification performance was not statistically significantly affected by inclusion/exclusion of demographic variables from the model. The Elastic net regression selected regions implicated in emotion and reward processing, as well as those involved in visual, auditory, and somatomotor processing. UD was associated with both reduced and increased cortical myelin in different areas, relative to HC. These results demonstrate that aberrant levels of cortical myelin may be a biomarker of unipolar depression.

Analyses examined 349 cortical parcels and 3 demographic variables, of which only 90 parcels and IQ were ever selected by the elastic net regression. The selection of IQ is unsurprising because the UD group had a higher mean IQ score than HC (Table 1). However, the model fit did not significantly change when no demographic variables were included. Of the 90 parcels selected at least once, 33 were selected more frequently than expected by chance based on permutation analyses.

There is abundant evidence of cognitive control and executive function deficits in depression⁵, and the anterior cingulate and orbitofrontal cortex play key roles in reward processing⁶⁴, which is also disrupted in depression¹¹.The 33 selected parcels included 13 parcels in regions important for executive functioning and cognitive control, located in the anterior cingulate (left a24pr, left 24dv, left p32, and left s32), orbitofrontal cortex (left 11l, left p10p, and left 47m), posterior cingulate (right RSC, and left d23ab), frontal operculum (right FOP1 and left FOP4), and parietal cortex (right POS2 and right 7Pm). UD was largely associated with reduced cortical myelin in these regions, except in one region of the anterior cingulate (left p32), and the two regions in the parietal cortex (right POS2 and right 7Pm), where it was associated with increased myelin. The present results converge with mounting evidence of disrupted myelination²⁸, thickness¹², connectivity⁵⁵, and activation¹¹ of these regions in those with depression.

In addition to regions traditionally reported in neuroimaging studies of depression, 19 of the identified parcels play roles in visual, somatomotor, and auditory processing. Visual processing regions included extrastriate regions (right V3A, left V4t, and left MT), ventral stream regions (right FFC and right TE2p), and dorsal stream regions (left MST, right PH, left 7PC, right 7AL). Auditory processing regions included regions in the auditory cortex (left LBelt, right LBelt, left RI, and right RI) and regions implicated in language (left STSdp, left FOP2, and right 52). Somatomotor regions included regions implicated in somatosensation (right 5m and right OP1) and a region in the premotor cortex (left 6a). UD was largely associated with reduced cortical myelin in auditory and somatomotor regions, except for two auditory regions (right 52 and right RI) where it was associated with increased myelin. In contrast, participants with UD had greater myelin in 6 visual regions and had reduced myelin in 3 regions (left V4t, left FFC, and left 7AL). While depression is not typically considered a disorder of dysfunctional sensory processing, somatic symptoms in depression have been well described⁶⁵, and there are reports of both visual⁶⁶ and auditory⁶⁷ processing deficits in depression. Our results contribute to a growing body of literature documenting associations of depression with altered structure¹² and functional connectivity^68,69 in sensory regions. These findings may point to disrupted information transfer in sensory regions a part of the underlying neurobiological etiology of depression.

Broadly, our results contribute to the consensus that the neurobiology of regions that play important roles in cognitive processes that are disrupted in depression may contribute to the etiology of the disorder. It is notable that in contrast to studies of the major white matter tracts, where depression is associated with lower integrity¹⁵, we observed both decreased and increased cortical myelin. This was particularly true in visual regions, where diagnosis was associated with increased cortical myelin in the majority of regions. Intriguingly, there are several examples in our results where proximal or similar regions had effects in opposing directions, including in the anterior cingulate (left p32 was increased, while left a24pr and left 24dv were decreased), the superior parietal cortex (left 7PC was increased and right 7AL was decreased), the retroinsular cortex (right RI was increased and left RI was decreased), and in visual processing regions (left MT and left MST were increased, and left V4t, which they border, was decreased). While the relationship between cortical myelin and functional activation or connectivity remains under-explored, these observations suggest that cortical myelin imbalance, rather than a global reduction, could drive some of the observed functional differences in depression, such as disrupted network integration^69,70 and patterns of both hypo- and hyper-connectivity^71,72.

Intriguingly, 84% of participants were consistently classified either quite accurately (>80% n=53) or inaccurately (<20% n=19) across nested cross-validation folds. In contrast, in permutation analyses where variables were retained, participants were classified at chance (50%). This suggests that inaccurate classification was largely not due to the algorithm guessing, but rather to the presence of brain features which reliably led to certain participants appearing like they were members of the other group. It was hypothesized that participants who were consistently inaccurately classified also had demographic or clinical features that were closer to members of the other group. For instance, perhaps UD participants who were inaccurately classified had fewer symptoms, or took different medications, from those who were accurately classified. Post-hoc analyses did not find any association between demographic or clinical variables and classification accuracy.

It is possible that inaccurate classification of HC participants as UD could reflect the presence of environmental or genetic risk factors that are not captured by clinical measures. As a preliminary exploration of this hypothesis, analyses tested the classification of a participant who was not used in primary analyses, as they experienced the onset of UD during the course of the study. This participant completed two MRI scans, at 12 months and 6 months prior to UD onset, and was classified as UD on both, despite not meeting criteria for a UD diagnosis at the time of scan. This preliminary result suggests that the pattern of cortical myelin in frontal, sensorimotor and extended visual cortices may be a biological marker of risk for UD diagnosis in the future. Further longitudinal studies are needed to test this hypothesis.

While the present results demonstrate that the pattern of cortical myelin is disrupted in UD, the cause of myelin disruption remains unknown and may include environmental, genetic, and other factors. However, post-hoc analyses suggest that myelin disruption is not due to demographic differences between the groups, nor due to clinical features of UD, such as medication use, lifetime depression severity, or illness duration. Disrupted cortical myelin may reflect other risk factors for UD. For instance, sleep disturbance, which is a well-established risk factor for UD⁷³, was recently shown to correlate with cortical myelin in several of the same regions found in the present study, including the cingulate, orbitofrontal cortex, and middle temporal cortex³². Another prominent hypothesis is that the source of white matter differences in depression is partially attributable to stress¹⁵. There is evidence that stress is associated with lower integrity of white matter tracts in humans^74,75 and in rodents and non-human primates it has been found to cause microstructural alterations to white matter^76,77. It is intriguing to note that remyelination can occur after stress as well⁷⁸, and there is evidence that it results in altered patterns of myelination across the cortex⁷⁹. This could potentially explain our observations of both decreased and increased cortical myelin in unipolar depression.

One limitation of this work includes the need to replicate our findings in an independent sample. To partially address this limitation, our analyses used robust machine learning methods involving model testing using held-out samples. Nested cross-validation was used to do model parameter selection without over-fitting, and cross-validation was used to estimate model accuracy, in which a model was iteratively trained and tested on non-overlapping participant subsets. Highly significant classification accuracy was achieved on held-out participants not used to train the model. The major strength of this approach is that it helps to reduce model bias, which occurs when the same participants are used to train and test a model⁸⁰. It has been suggested that model performance with nested cross-validation is close to the accuracy that would be achieved on fully independent data⁸¹. It is also notable that the observed model accuracy is greater than the association of any of the individual parcels selected. The second limitation concerns the reliability of the T1w/T2w metric. While it is reported to be relatively insensitive to noise⁴⁹, there have not been any studies of the test-retest reliability of estimates of cortical myelin. Future research will be needed to verify that the T1w/T2w ratio has adequate reliability. The third limitation concerns measurement noise due to susceptibility artifacts. Several regions with documented relevance to depression, including the bilateral hippocampus, entorhinal cortex, and posterior orbitofrontal cortex complex, were not included in the present analyses, as these regions showed an excess of between-person variability (see Methods). Future research should explore the association of cortical myelin in these regions with unipolar depression.

In summary, cortical myelin can distinguish participants with UD from HC, even when clinical and demographic variables are not included in analyses. Regions that were most important for this classification include several that play key roles in reward and emotion processing as well as a host of regions important for sensory processing. This result highlights that the association of UD with sensory processing bears further investigation. Notably, UD was associated with both decreased and increased cortical myelin, suggesting that observations of reduced integrity of major white matter tracts in UD may not fully extend to the cortex. These results suggest that cortical myelin holds promise as a biomarker of unipolar depression and may be an early predictor of risk for this disorder.

Funding Acknowledgements

This work was supported by a grant from the National Institute of Health R01MH114870 to A.M., and Y.O.H was supported by P41EB019936 to the Center for Reproducible Neuroimaging Computation (PI: Kennedy). D.A.A.B, Y.O.H., S.S., R.R, S.I., and A.M.: declare no conflict of interest. H.A.S: receives royalties from Wolters Kluwer, royalties and an editorial stipend from APA Press, and honorarium from Novus Medical Education.