Major depressive disorder on a neuromorphic continuum

Introduction

Major depressive disorder (MDD) is a devastating worldwide mental disorder1. Patients with MDD are clinically diagnosed according to the DSM-5-TR when a patient exhibits at least five designated polythetic symptoms2. Patients with remarkable heterogeneous clinical manifestations may have co-occurring subtypes of symptoms, and therefore, understanding the biotypes of shared brain signatures may improve our understanding of etiological mechanisms3. However, such a categorical model does not consider alternative latent effects4. Indeed, depression has been more recently conceptualized as an unbroken “spectrum” (also referred to as factors), rather than as a “category”5. This conceptualization raises the possibility that neuroimaging can be used to identify more granular disease “patterns”, accounting for individual and specific system reorganizations, which may be more readily translatable to individualized patient care.

However, highly nosological and neurobiological heterogeneities exist among MDD patients, including subjective symptoms6, neuroimaging intermediate phenotypes7,8, and genotypes9. Analyses of interindividual differences in MDD populations have therefore been conducted by numerous investigators, as shown by recent attempts to stratify individuals with MDD based on brain neuroimaging or cognitive skills10. Prior neuroimaging-based categorical subtyping studies in MDD have assumed that each participant belonged to a distinct categorical subtype based on a shared neuroimaging features base11,12,13,14,15,16. Progress toward translating knowledge into clinical practice has remained challenging17. One possible explanation may involve the clustering analytical approaches used to separate categories18. Indeed, categorical models merely highlight common brain features, whereas they disregard continuous individual variations within subtypes19. In contrast to a “winner-take-all” assumption, multiple domain studies suggest that MDD is a continuous entity instead of a collection of categories, each with clearly defined boundaries20,21,22.

Continuous variations across individuals with MDD have been observed with varying degrees across various symptom domains (exo-phenotypes)23,24. In parallel, evidence from genetic variants has suggested that further insights into depression genotypes can be obtained when this within-group variation is considered25,26,27. Both exo-phenotype and genotype findings have confirmed the continuous variations across individuals, leading to high variance, even within subtypes28,29,30. Thus, we hypothesize that it is more plausible that individuals with MDD exhibit distinct varying degrees of expressions of endo-phenotypes. An alternative approach to conceptualizing MDD heterogeneity based on brain neuroimaging-derived phenotypes is dimensional methods, such as the data-driven Bayesian framework, to assess the co-expression patterns of observable and latent abnormalities within each individual that may reflect co-existing pathologies. The mathematical framework, latent Dirichlet allocation (LDA), facilitates the possibility that multiple latent factors are expressed to varying degrees within a patient31. LDA has been successfully used to decompose whole-brain structural or functional features into abnormal patterns in various brain disorders32,33,34,35,36,37. For example, a patient’s abnormal pattern may be 80% for factor 1 and 20% for factor 2, while another patient’s abnormal pattern may be 60% due to factor 1 and 40% due to factor 2. Given that multiple contributors are not mutually exclusive and, therefore may influence heterogeneity in MDD38, we believe that it is more biologically plausible that an individual expresses multiple disease factors (i.e., one or more) to various degrees (continuum), rather than by assigning each individual to a single subgroup. This approach is particularly relevant to the hypothesis in the current study.

Magnetic resonance imaging (MRI) is a beneficial technique because it is a reliable, noninvasive measurement, which can be easily used for diverse individuals with brain disorders to quantify brain morphometries, such as cortical thicknesses and volumes39,40. MDD is associated with widespread brain structural changes41. Cortical thickness is thought to reflect the density and arrangement of neurons, neuroglia, and nerve fibers42. Furthermore, meta-analyses have reported cortical thickness and subcortical volumetric abnormalities in MDD8,43. Recent studies have reconciled previously inconsistent findings of either decreased or increased morphometric characteristics by showing that both patterns co-exist, with each affecting distinct anatomical regions7,8,43,44. However, these studies involve the “one-size-fits-all” theory, which is inconsistent with the observations of clinicians. Thus, when considering interindividual variations across MDD patients, we sought to provide a detailed characterization of the nature and spatial extent of brain MRI-derived phenotypes in these patients.

This study investigated latent disease factors (spatial pattern) and continuum factor expressions in individuals with MDD. To this end, we used a Bayesian model to decompose structural MRI features (whole-brain cortical thickness and subcortical volume) in a multisite cross-sectional cohort (n = 2294). We showed that (1) it is reasonable to identify expressions of multiple abnormal patterns rather than assigning each individual into a single subgroup of MDD patients; (2) individual factor expressions of patients with MDD were stable over time in a longitudinal cohort and could be replicated in an independent replication cohort; (3) individual factor expressions of MDD patients showed disease specificity when compared with three other mental disorder controls, in an independent transdiagnostic cohort; and (4) these latent disease factors (spatial pattern) robustly predicted various treatment responses in the longitudinal cohort. Latent factors in MDD, therefore, provided the means to both dissect heterogeneous categories and to better understand the psychopathology of MDD patients.

Results

Study design

To identify reliable latent factors of MDD, we included MDD patients (n = 1079) and healthy controls (HCs, n = 1215) from six sites in discovery and replication cohorts. The demographic characteristics of each scan site, after various exclusion criteria based on data quality (Methods) are presented in Supplementary Table 1. One site [i.e., First Hospital of ShanXi Medical University (FHSXMU)] included both cross-sectional and longitudinal samples. We thus used this site data as a longitudinal cohort to avoid repeated usage of treatment responses. The study design is shown in Fig. 1. Using factor specificity analysis, we enrolled 380 cases with three other psychiatric disorders from a transdiagnostic cohort: patients with bipolar disorder (BD, n = 103), schizophrenia spectrum disorder (SSD, n = 147), and generalized anxiety disorder (GAD, n = 130) (Supplementary Table 2). We examined factor robustness and assessed the performance of latent disease factors to predict treatment responses based on five longitudinal datasets (n = 122), comprising individuals with repetitive transcranial magnitude stimulation (rTMS, n = 16), electroconvulsive therapy (ECT, n = 49), neurofeedback (n = 9), fluoxetine (n = 14), and ketamine (n = 34) (Supplementary Table 3). The distribution of ages is presented in Supplementary Fig. 1.

Fig. 1: Study design.
figure 1

a This study aimed to identify latent factors of patients with major depressive disorder (MDD) based on cortical thickness (CT) and subcortical volume. We decomposed patients with MDD into three latent disease factors (spatial pattern) and obtained factor compositions (individual expression) based on multisite cohorts. We then replicated the disease factors (spatial pattern) in an independent cohort. b To test factor specificity, we enrolled other three psychiatric disorders, including bipolar disorder (BD), generalized anxiety disorder (GAD), and schizophrenia spectrum disorder (SSD). c We also used longitudinal data to investigate whether the factor compositions (individual expression) were stable over time. d Based on longitudinal data, we also determined the clinical efficacies of the latent factors. Cartoons were created in BioRender. Liao, W. (2025) https://BioRender.com/j18v472.

Full size image

Latent abnormal factors in MDD patients

Considering inherent continuous variation among individuals with MDD, we used the dimensional Bayesian LDA model, an unsupervised machine learning technique originally devised for text mining, to decompose the structural MRI of patients with MDD. This decomposition framework combines categorical and dimensional models, enabling us to test whether MDD patients expressed one or more latent abnormal factors. After controlling for data quality, we used structural MRI (M features) of MDD patients (i.e., defined as a discovery cohort, n = 928) (Supplementary Fig. 2). After inputting n×M into the model, where n was the number of patients and M was the number of cortical vertices and subcortical structures, we obtained the probabilistic abnormal map of the particular disease factor (i.e., disease load, and spatial pattern) and the factor composition of an individual (i.e., factor probability for each individual). To identify the latent factors, K, we evaluated a range of latent factors of K from 2 to 5. Because the 4-factor and 5-factor models were unstable, we did not experiment with more factors (Supplementary Fig. 3a). The 3-factor model achieved the largest spatial separation between latent factors (Supplementary Fig. 3b), so we next focused on the results of the 3-factor solution.

Three factor-specific patterns across cortical and subcortical cortices are shown in Fig. 2. After performing false-discovery rate (FDR) corrections, we found that Factor 1 was associated with MDD-related increased cortical thickness in sensory cortices. In contrast, MDD-related decreased cortical thickness was located in the anterior insula and orbitofrontal cortex. Factor 2 was dominated by patterns of decreased cortical thickness in the cingulo-opercular network (CON)45. Factor 2 also showed a decreased subcortical volume, particularly in the amygdala, nucleus accumbens, and hippocampus. Furthermore, Factor 3 was mainly characterized by patterns of increased cortical thickness in social and affective brain systems46,47.

Fig. 2: Patterns of three latent disease factors in major depressive disorder (MDD).
figure 2

ac represent probabilistic maps of Factor 1, Factor 2, and Factor 3, respectively, in MDD patients. The brighter color indicates a higher probability at the cortical vertices and subcortical structures for a given factor [i.e., Pr(Pattern|Factor)]. The presented values are arbitrary units.

Full size image

Examination of abnormal factor compositions among MDD patients revealed that most patients expressed multiple latent factors, rather than a single factor (Fig. 3a), indicating a continuum of individual variabilities. In addition, no single site showed predominantly a single factor, suggesting that site differences did not drive latent factors. To validate the 3-factor patterns derived from the discovery cohort, we used the Bayesian LDA model for the validation cohort (i.e., FHSXMU dataset), and found that the 3-factor model exhibited similar patterns across the whole brain between the discovery and validation cohorts (Factor 1: R = 0.58, P < 0.0001; Factor 2: R = 0.50, P < 0.0001; Factor 3: R = 0.52, P < 0.0001; Supplementary Fig. 4).

Fig. 3: Factor composition and specificity.
figure 3

a Factor compositions in major depressive disorder (MDD) patients. Each patient is represented as a dot in the triangle, with its barycentric coordinate representing the factor composition expressed as posterior probability [Pr(Factor|Patient)]. A dot located close to the corners predominantly expresses a given factor (i.e., dots with red, green, or purple colors), whereas those located toward the centroid of the triangle express the combination of three factors. b Specificity of factors. The severity of factor expression [z-score weighted by disease factors Pr(Pattern|Factor)] was projected into a three-dimensional space and colored by groups, including MDD patients (n = 928), healthy controls (HCs, n = 1104), and disease controls composed of patients with bipolar disorder (BD, n = 103), generalized anxiety disorder (GAD, n = 130), and schizophrenia spectrum disorder (SSD, n = 147). Data are presented as mean values and a confidence interval of 75%.

Full size image

After validating the robustness of factor patterns and composition, we also evaluated factor specificity using the other three mental disorders. For each individual, we calculated the average z-score weighted by the posterior probability for each disease factor obtained from MDD patients. We found that MDD factors were not expressed in HCs, and were expressed in opposite order in the other three psychiatric disorders, including BD, GAD, and SSD (Fig. 3b). In addition, we identified latent factors in the three disorders to determine whether abnormal patterns were implicated in each of the factors with discovery sample (Supplementary Fig. S5), and measured the percentage of overlapping voxels across factors between discovery and the three cohorts (Supplementary Table 4).

Associations between abnormal factors and neurotransmitter receptors

Given that neurotransmitter dysregulation plays an important role in depression, we next sought to quantify the associations between neurotransmitter dysregulation and the abnormal factors used in the study. Neurotransmitter receptors drive synaptic plasticity, modify neural states, and ultimately form network-wide communications. We included 19 distinct neurotransmitter receptors and transporters across nine different neurotransmitter systems, which were collated from open positron emission tomography (PET) tracer images48. After assigning the receptor distribution and three disease factors to cortical DK-68 and subcortical brain regions (Fig. 4), we used a multiple linear regression model to combine all 19 neurotransmitter receptors and transporters using the relaimpo package in R. The model explained more than 85% of the variances for all three abnormal factors in MDD patients (Factor 1: Adjusted R2 = 0.88, F(19, 62) = 32.59, P = 2.2e-43, Fig. 4a; Factor 2: Adjusted R2 = 0.92, F(19, 62) = 49.29, P = 4.8e-50, Fig. 4b; Factor 3: Adjusted R2 = 0.89, F(19, 62) = 35.03, P = 1.5e-44; Fig. 4c). However, the results differed from the significant individual neurotransmitter receptors and transporters among the three abnormal factors. We found that the major predicted neurotransmitter receptors and transporters in Factor 1 were norepinephrine (P = 0.0002, FDR-correction), cannabinoid (P = 0.0002, FDR-correction), and 5-HT2A (P = 0.0002, FDR-correction). The norepinephrine (P = 0.00003, FDR-correction), 5-HTT (P = 0.0002, FDR-correction), and mGluR5 (P = 0.0006, FDR-correction) were primarily associated with Factor 2 distribution. For Factor 3, the serotonin system, including 5-HT2A (P = 0.0006, FDR-correction), 5-HT4 (P = 0.02, FDR-correction), and 5-HTT (P = 0.0006, FDR-correction), mainly explained the cortical-subcortical distributions. Together, these findings revealed that the three factors corresponded to a diverse and distinct set of neurotransmitter receptors and transporters, which provided data regarding heterogeneity in MDD patients while identifying potential mechanisms involving distinct abnormal factors.

Fig. 4: Contributions of neurotransmitter receptors/transporters to latent disease factors.
figure 4

ac were the observed and predicted disease factors (spatial maps) 1, 2, 3, respectively, which were projected on the DK-68 atlas and subcortex. The right panel plots the predicted versus observed values (n = 82; a: P-values = 2.2e-43; b: P-values = 2.2e-43; c: P-values = 1.5e-44). The statistical significance was obtained by the F-test. d The significant neurotransmitter receptors or transporters in the predicted model for each factor pattern (PFDR < 0.05) are shown.

Full size image

Individual predictions of clinical outcomes

We next aimed to determine whether the expressions of abnormal factors remained stable over time, using longitudinal MRI data of MDD patients across four centers (n = 122). We used probabilistic abnormal maps estimated from MDD patients in the discovery cohort to infer factor compositions of patients from longitudinal data (i.e., pretreatment and posttreatment MRI data) using the standard variational expectation-maximization (VEM) algorithm. To determine whether abnormal factors reflected different disease stages rather than abnormal subtypes, we compared factor compositions of MDD patients across longitudinal data. We found that MDD patients before and after treatment exhibited more than one latent factor (Fig. 5a). The factor probabilities were positivity correlated and highly consistent (Factor 1: R = 0.75, P < 0.0001; Factor 2: R = 0.92, P < 0.0001; Factor 3: R = 0.89, P < 0.0001; Fig. 5b), suggesting that the factor compositions of MDD patients were stable, despite disease progression and clinical treatments.

Fig. 5: Stability of factor compositions across the longitudinal cohort.
figure 5

a We inferred the factor compositions for each individual in the longitudinal cohort, based on the constructed Bayesian model from the discovery cohort. Each dot represented an individual in a longitudinal cohort. Each color represents a clinical treatment, and five treatment options were included in this study, including electroconvulsive therapy (ECT), ketamine, repetitive transcranial magnetic stimulation (rTMS), fluoxetine, and cognitive behavioral therapy (CBT). Left and right panels were the factor compositions of pretreatment and posttreatment, respectively. b Stability of factor compositions over time. Each dot represents a participant. For each plot, the x and y-axes represent the probabilities of a factor at baseline and after baseline, respectively. P-values are obtained by Pearson’s correlation analysis with two-sided (n = 122; F1: P = 6.4e-23; F2: P = 6.7e-51; F3: P = 1.5e-43). Different colors denote different treatment options.

Full size image

We finally determined whether the findings of abnormal factors were related to clinical efficacy. We tested whether disease severity predicted depression improvement in physical interventions (rTMS, ECT, and neurofeedback) and the pharmacological treatments (fluoxetine and ketamine) group by using the partial least squares regression (PLSR) method. There was a significant positive relationship between the observed and predicted improvements in depression severities (physical interventions: R = 0.26, Pperm = 0.01; Fig. 6a, left panel; pharmacological treatments: R = 0.33, Pperm = 0.04; Fig. 6b, left panel), suggesting a potential clinical application of using the identified factors. To further investigate the relationships between clinical outcomes and abnormal factors, we grouped the MDD patients into three scales based on treatment outcomes (i.e., non-responders, partial responders, and responders). There was no statistical difference in demographic information among the three groups in physical interventions (age: F = 1.67, P = 0.20; sex: χ2 = 0.06, P = 0.97), and in pharmacological treatments (age: F = 0.35, P = 0.71; sex: χ2 = 0.20, P = 0.90). We next assigned MDD patients to the three factors according to their factor compositions, using a “winner-take-all” approach. In physical interventions, 24% of MDD patients were non-responders, 24% were responders, and 54% were partial responders (Fig. 6a, upper panel). In pharmacological treatments, 23% of MDD patients were non-responders, 38% were responders, and 40% were partial responders (Fig. 6b, upper panel). In addition, non-responders showed high expressions of Factor 1 in the physical therapies and pharmacological treatments group (Fig. 6a, b, lower panel). Partial responders and responders showed high expressions of Factor 2 in physical therapies, and Factor 3 in pharmacological treatments (Fig. 6a, b, lower panel). However, there were no significant statistical comparisons.

Fig. 6: Clinical outcome predictions based on latent factor patterns.
figure 6

a Clinical improvement prediction of physical intervention. b Clinical improvement prediction of pharmacological treatment. Left panel: A plot of observed and predicted clinical improvements based on the three latent factors. Each dot denotes an individual in the longitudinal cohorts (a: n = 74; b: n = 48). The significance (Pperm) of Pearson’s correlation coefficients are evaluated using the permutation test (one-sided). Upper panel: the correspondence between response degrees for clinical treatments and factor patterns. Green denotes non-responders, red denotes partial responders, and purple denotes responders. Lower panel: in section a, the distribution of factor compositions for non-responders (F1: n = 1; F2: n = 13; F3: n = 3), partial responders (F1: n = 1; F2: n = 11; F3: n = 5), and responders (F1: n = 3; F2: n = 36; F3: n = 1); in section b, the distribution of factor compositions for non-responders (F1: n = 3; F2: n = 6; F3: n = 2), partial responders (F1: n = 1; F2: n = 13; F3: n = 4), and responders (F1: n = 2; F2: n = 12; F3: n = 5). Data are presented as mean values +/− standard deviations.

Full size image

Discussion

This study identified three latent abnormal factors in MDD patients, using a Bayesian model of cortical thickness and subcortical volume maps. While not using the conventional case-control analyses and prior subtyping models that assigned each individual to a single subtype, the LDA approach addressed the continuum of interindividual variabilities in MDD patients, by estimating the factor composition of multiple abnormal factors for each patient. Specifically, each factor was expressed to different degrees across MDD patients, and each was associated with distinct neurotransmitter receptors or transporters collated by using open PET sources from another 1200 healthy individuals. By identifying unique relationships between the main structural biomarkers of MDD patients, which were disease-specific compared with three other psychiatric disorders and stable over time across the longitudinal sample, this integrative analysis increased our understanding of the complex landscape of this condition. Notably, the three latent factors contributed to predicting clinical outcomes, further illustrating the potential ability of dimensional modeling to improve clinical efficacies.

The increasing popularity of analyzing interindividual variations in MDD populations has been demonstrated by recent efforts to categorize MDD patients according to the degree of clinical symptoms and brain neuroimaging abnormalities41. However, many symptom domains manifest as continuous variations across individuals26, even within subtypes49. Given that interindividual variability in MDD probably reflects varying degrees of factor expressions and related mechanisms, we used the polar LDA technique, which was chosen in consideration of the increasingly acknowledged fact that variations in neuroimaging signatures in psychiatry constitute a dimensional continuum, to disentangle MDD-specific neuroanatomical variation from variation shared with HCs. Our 3-factor model revealed the following: a complex pattern of increased unimodal cortical thickness and decreased orbitofrontal and insula cortical thickness; a cortical pattern associated with decreased cortical thickness in the CON and a subcortical pattern in the amygdala, nucleus accumbens, and hippocampus; and an increased cortical thickness factor associated with the social and affective cortices. Our analyses, integrating neuroimaging and molecular factors, therefore revealed neurotransmitter receptors/transporters most correlated with 3-factor patterns, which were implicated in the pathophysiology of MDD patients.

Consistent with our previous findings50, some MDD patients exhibited increased cortical thicknesses within unimodal cortices, in Factor 1. Using the data-driven analysis across the whole brain, Anderson et al. also reported the surprising pattern of increased cortical thickness in visual cortices in MDD patients, relative to controls, which could be missed by hypothesis-driven examinations of selected brain regions51. In addition, they found a similar neuroimaging pattern to that of this study with Factor 1, where MDD patients exhibited reduced global brain connectivity in the orbitofrontal cortex, but increased global brain connectivity in the visual cortex. More importantly, increasing studies have reported that abnormal functions or structures of visual cortices were related to depression and antidepressant efficacies52. Other studies reported two converse MDD neuroimaging patterns related to default mode, limbic, sensorimotor, occipital, and subcortical regions12,16,53,54. However, in contrast to these findings that patients with MDD were clustered into a specific subtype based on a shared neuroimaging pattern, most patients expressed multiple factors, suggesting a mosaic of subtypes within individuals with MDD. In the present study, we hypothesized that abnormalities in MDD patients might be associated with decreased cortical levels of the inhibitory neurotransmitter GABA, in the occipital cortex55,56. Specifically, our multivariate analyses found that the Factor 1 pattern was associated with the GABAA receptor distribution map, crucial for maintaining neural balance57. In addition, the cannabinoid receptor 1 map was highly correlated with the Factor 1 pattern and was expressed in structures implicated in depression, such as the frontal and limbic brain regions58. Recent studies have also suggested that cannabinoid receptor 1 was a promising therapeutic target for anti-depression and anti-anxiety treatments59. The key brain regions in Factor 1, including orbitofrontal and occipital cortices, may become novel target areas for treating depressed patients52,60.

In the present study, Factor 2 was characterized by decreased cortical thickness within CON and subcortical volume in the amygdala, hippocampus, and nucleus accumbens (NAc). At the neuronal level, the CON closely resembled the network most strongly associated with behavioral and emotional regulation61. Specifically, the CON was structurally and functionally altered in MDD patients62. The dorsal anterior cingulate cortex is considered a core hub within the CON. Previous studies have reported that the CON was a key area involved in the pathology of MDD, which was found using functional MRI, metabolic studies, and perfusion approaches to be altered in MDD patients63. In addition, the NAc is a critical region for reward and motivation, and its dysregulated activity is a hallmark of MDD. Notably, a previous study reported that deep brain stimulation to the NAc improved depressive-like behaviors in an MDD mouse model64, suggesting the critical target role of NAc in MDD patients. In the present study, the norepinephrine transporter was the neurotransmitter making the most contribution to Factor 2. It has long been known that norepinephrine plays a crucial role in brain processes that are responsible for the development of symptoms associated with severe depressive illnesses65. Furthermore, norepinephrine transporters, binding in the locus coeruleus, reflect the dynamics of norepinephrine and are decreased in postmortem tissues from patients with MDD66, suggesting that this neurotransmitter may be a target for tricyclic antidepressants.

In the present study, in contrast to Factor 2, Factor 3 mainly exhibited increased cortical thicknesses within the social and affective-associated brain regions. The temporoparietal junction (TPJ) is essential for multiple cognitive and social functions and may act as a critical node for integrating internal and external information. Numerous neuroimaging studies have reported TPJ abnormalities in MDD patients, compared to healthy controls67,68. Furthermore, a Factor 3 pattern has been correlated with specific neurotransmitter systems, such as serotonin receptors and dopamine transporters. The serotonin hypothesis of depression, whose abnormalities result in depression, has been studied for decades, providing essential justification for the use of some antidepressants. Although the serotonin hypothesis is controversial69, molecular imaging has revealed that the system is disturbed, and acute tryptophan depletion and reduced plasma tryptophan during depression suggest a function for 5-HT in patients who are susceptible or experiencing depression70. In this study, we hypothesized that Factor 3 might influence clinical aspects through the serotonin system.

With regards to the prediction of clinical outcomes, previous studies have linked subnetworks, such as the triple and sensory networks, to depression symptom relapse after clinical treatments3. In the present study, integrating the 3-factor patterns predicted improvements in depression severity in MDD patients, showing that imaging-based disease severity (i.e., disease load weighted factor composition) was related to state change. In this model, non-responders showed the highest factor compositions within Factor 1, when compared with partial-responders and responders, which suggested the adverse neuroplasticity effect of complex cortical abnormal patterns11. Furthermore, analysis of longitudinal MRI scans revealed that factor compositions were stable over time, indicating that individuals were not progressing from one factor to another, and their factor compositions were unaffected by their mood state. The observed heterogeneity in MDD might correspond to disease expressions, rather than disease stages34.

Our study had several limitations. First, because of its multi-center design, we could not associate the three latent factors with clinical symptoms. The MDD patients had various cognitive impairments that were not assessed in the current retrospective study. Further analyses should unify the assessments of clinical symptoms with combined detailed cognitive performances, to better understand the complex relationships between the neurophysiological basis and clinical presentations of MDD patients. In addition, medication status is also a confounding factor in neuroimaging analyses. However, due to the lack of information, we did not explore the medication effects on our results. Future studies should, therefore, conduct a detailed analysis of clinical variables and possibly investigate the relationships between clinical variables and neuroimaging analyses. Second, the sample size of the longitudinal cohort was too small across clinical treatment types. In addition, the period of treatment strategies varied. The period and treatment strategy differences may influence the predictions of individual treatments. We thus separated the patients into two groups: physical interventions and pharmacological treatments, to predict individual treatment. However, future studies should include more patients with clinical treatment in each treatment strategy, to better establish prediction models for non-responders, partial responders, and responders. Third, previous studies have reported that in-scanner head motion is an essential noise source in cortical thickness and gray-matter volume studies71,72,73,74. Given the lack of an effective method to quantify head motion from 3D-T1w images, we used Freesurfer’s Euler number as covariates in our model. Euler number is considered to be consistently correlated with manual ratings across samples75 and was reduced in motion-affected group76. Future studies should also use volumetric navigators in the T1 sequence to measure head motion. Finally, information on neurotransmitter receptors/transporters has been mainly obtained from patients’ brains in Europe, and may not entirely reflect the physiology of the Asian participants in our study. Future studies need to integrate multi-omics data from the same samples to verify these findings.

In summary, using a polar LDA approach, our study identified three latent abnormal factors in MDD with increased and decreased cortical thicknesses or subcortical volumes. The factors were associated with distinct neurotransmitter systems. Furthermore, our findings showed that each patient expressed multiple abnormal latent factors to various degrees, rather than a single factor, which could improve the prediction of depression severity. Therefore, each patient’s unique factor composition might be helpful for future biomarker development and clinical treatments. Together, our findings revealed that individuals with MDD were organized along continuous dimensions that affect distinct sets of brain regions, which might be useful for the future identifications of biomarkers.

Methods

Participants

MRI data used in this study were obtained from nine research groups (Supplementary Methods). In the discovery cohort, we included 928 patients with MDD (533 females and 395 males) and 1104 healthy controls (652 females and 452 males). Institutional review boards approved the study procedures for the participating institutions, and written informed consent was obtained from all participants.

Structural image processing

All structural MRI data of participants were obtained using 3 T MRI scanners. The detailed data acquisitions for each site are described in Supplementary Table 5. The T1w images were automatically preprocessed using the FreeSurfer version 6.0 (https://surfer.nmr.mgh.harvard.edu/fswiki). Briefly, following intensity normalization, non-brain tissues were removed from individual T1 images. The gray and white matter were segmented and then used for cortical reconstruction. The white surface (defined as the gray/white boundary) and pial surface (defined as the gray/cerebrospinal fluid boundary) were then identified. Cortical thickness was calculated as the closest distance from the gray/white boundary to the gray/CSF boundary at each vertex, resulting in a cortical thickness map for each participant. The individual maps of cortical thicknesses were then warped and registered to standard space (i.e., 4661 vertices in fsaverage4), thus, enabling matching of cortical locations among individuals across the whole surface. In addition, the subcortical volumes of the bilateral thalamus, caudate, NAc, pallidum, putamen, amygdala, and hippocampus were also obtained. The cortical thicknesses at each vertex and subcortical volume were combined (4675 features in total) and used for the LDA analyses.

Quality control of MRI data

The quality of T1w images was assessed using the MRIQC. T1w data were excluded if the scan was marked as an outlier (1.5 × the inter-quartile range in the adverse direction) in the following six image quality metrics: coefficient of joint variation, contrast-to-noise ratio, signal-to-noise ratio, foreground-to-background energy ratio, entropy focus criterion, and full-width half maximum smoothness. Detailed information about these image quality metrics is found in the MRIQC’s documentation.

Normative modeling

After controlling image quality, we used normative modeling to obtain individual-specific cortical vertices’ thickness and subcortical volume deviation maps concerning an underlying model of normative expectations for each vertex variation. Normative modeling was run using the Predictive Clinical Neuroscience toolkit package (PCNtoolkit, version 0.26, https://pcntoolkit.readthedocs.io/en/latest/). To increase the sample size in constructing a normative model, we included healthy controls from Human Connectome Projects (n = 1113; age: 22–37 years) and all five sites. However, multisite neuroimaging data may introduce artefactual variance, impacting subsequent analyses. Thus, to account for these site-related effects, we used the recommended warped Bayesian Linear Regression (BLR) implemented in PCNtoolkit77, based on Bayesian linear regression with likelihood warping, to predict cortical thickness and subcortical volumes, including age, sex, site, and Euler number (summarizes the topological complexity of the reconstructed cortical surface and has been shown to reliably detect reduced image quality in 3D-T1w scans75) as covariates78.

Briefly, for each brain region of interest (cortical thickness or subcortical volume in this study), the predicted (hat{y}) was assumed to be a result of a linear combination of the B-spline basis function79:

$$hat{y}={w}^{T}{{varnothing }}left(xright)+{epsilon }_{s},ldots ldots ldots ldots ldots ldots .$$
(1)

where ({w}^{T}) is the estimated vector of weights, (varnothing left(xright)) is the basis expansion of the covariate vector x, including age, sex, site, and Euler number in this study, ({epsilon }_{s}=eta (0,{beta }_{s}^{-1})) is a Gaussian noise distribution for site s.

Deviation scores (i.e., z-scores) were calculated for patient i, and at each cortical vertex or subcortical region j, as:

$${z}_{{ij}}=frac{{y}_{{ij}}-{hat{y}}_{{ij}}}{sqrt{{sigma }_{{ij}}^{2}+{sigma }_{{nj}}^{2}}},ldots ldots ldots ldots ldots ldots .$$
(2)

where ({y}_{{ij}}) and ({hat{y}}_{{ij}}) are the actual and predicted cortical thickness and subcortical volume, respectively. ({sigma }_{{ij}}^{2}) is the estimated noise variance, reflecting uncertainty in the data, and ({sigma }_{{nj}}^{2}) is the variance attributed to modeling uncertainty. This model quantified the extent to which an individual’s cortical thickness or subcortical volume estimate deviated from the model prediction, given the uncertainty of the model.

Factor specificity analysis

To assess the specificity of latent factors in MDD patients, we quantified the expression of factors (i.e., disease severity) in HCs, BD, GAD, and SSD patients. For each individual, we defined the disease severity as the averaged absolute z-score across the top 5% from 4675 features for each disease factor weighted by its posterior probability [Pr(Pattern|FactorK), K = 1,2, …, Kmax]32. Before disease severity computation, the z-score was subtracted from the mean values and divided by the standard deviation to minimize site and disease effects.

Latent factors estimation

After obtaining the z-score data matrix (n patients × M brain features = 928 × 4675) in the discovery cohort, we identified latent MDD patient factors based on a Bayesian model named LDA31. The model used in the present study has been used to characterize latent disease factors in Alzheimer’s disease34, epilepsy32, and autism study33. The model assumed that each participant expresses one or more latent factor(s) related to distinct patterns of structural alterations. The z-scores were then converted to counts by multiplying them by 10 and rounding to the nearest integer. This was because data discretization might lead to some loss of information, so previous studies have suggested that a sufficiently large multiplicative factor (i.e., multiplied by 10) results in no loss of information32,33,34. Positive/negative counts denoted increased/decreased cortical thickness or subcortical volume in MDD patients relative to controls. Given the z-score data matrix and a specific number of factors, K, we obtained the factor composition, [Pr(Factor|Patient)] of each participant, indicating the continuum probability that a patient expressed a latent factor and latent disease factor, [Pr(Patterns|Factor)], indicating the spatial patterns of alterations of each factor.

We estimated the number of factors, K, from 2 to 5. The results showed robust Factors 2 and 3, and unstable Factors 4 and 5. Thus, we did not consider a larger number of factors. To determine which number, K, of latent factors (2 to 5 factors) was the best, the maximal spatial separation criterion was used80, which meant that the spatial distributions of the optimal K latent factors were maximally separated or overlapped minimally. The maximal spatial separation criterion was conducted as follows. First, Pearson’s correlation coefficients of all pairs of K latent factors were calculated. The separation values were then defined as the maximum correlation coefficient, and computed as follows:

$$1-{max }_{forall i,jin {1,ldots,K}}({corr}left(Pr left({Pattern} | {{Factor}}_{i}right)right),{corr}({Pattern|}{{Factor}}_{j})),ldots ldots ldots ldots$$
(3)

This procedure was performed for each K from 2 to 5. The 3-factor estimation had the highest separation value.

The bootstrapping procedure generated 100 samples from the patients’ data to assess confidence intervals of each factor-specific pattern of structural alterations, as suggested by a previous study32. The z-scores were then computed by dividing the pattern of structural alterations [Pr(Pattern|Factor)] by the standard deviation of 100 samples and converted to P-values. Multiple comparisons were corrected using an FDR of 0.01. The top 10% significant vertices were displayed to visualize the factor patterns better.

To assess the robustness of the latent factors estimated in the discovery cohort, we replicated the latent factor estimation in the replication cohort. The procedures for latent factor estimation were the same as those previously mentioned. The spatial correlation coefficient was then used to determine the similarity of structural alteration patterns of each latent factor, between discovery and replication cohorts.

Molecular fingerprints underpinning MDD latent factors

The brain-wide volumetric neurotransmitter receptor densities were estimated using 19 open PET-derived tracer images from a combined total of more than 1200 healthy participants across multiple studies. PET images were all normalized to the MNI-ICBM 152 non-linear 2009 and provided by Hansen et al. 48 and Markello et al. 81. These 19 receptors/transporters combined into nine neurotransmitter systems as follows82: dopamine (D1, D2, DAT), norepinephrine (NET), serotonin (5-HT1A, 5-HT1B, 5-HT2, 5-HT4, 5-HT6, 5-HTT), acetylcholine (α4β2, M1, VAChT), glutamate (NMDA, mGluR5), GABA (GABAA), histamine (H3), opioid (MOR), and cannabinoid (CB1). To decrease dimensions and noise effect in multiple correlation analysis, we rewarped the cortical DK-68 parcellation and 14 subcortical nuclei onto the volumetric PET image space. Next, neurotransmitter receptor profiles were parcellated to 82 brain regions and were z-scored, resulting in an 82 regions × 19 receptors matrix of relative densities.

A multivariate linear regression model combining 19 neurotransmitter receptors was used to identify molecular contributions to the three MDD latent factors83,84,85. The latent disease factors, [Pr(Pattern|Factor)], as responder and receptor matrix, as predictors were included in the relaimpo package (relative importance of regressors in linear models, version 2.2-5) in R. Relative importance metrics were used to address linear regression with multiple collinear regressors, and statistical significance was set at an FDR correction with P < 0.05.

Prediction of depressive symptom improvements in longitudinal cohorts

For longitudinal cohorts, the latent factor model that was first trained using the discovery cohort was used to determine the factor composition of both pretreatments and posttreatment in the longitudinal sample using the VEM algorithm33,34. The similarities of factor compositions between pretreatments and posttreatments were assessed by spatial correlation across all longitudinal MDD patients, to check if the factor composition was stable during disease progression. We further assessed the performance of latent factors to predict depressive symptom improvements in MDD patients. The PLSR is a regression technique that can behavior from brain features86. Leave-one-out cross-validation method was developed for this prediction model. Among the individuals, an individual was selected as a testing set, and the remaining individuals were used as a training set. To train the prediction model, we inputted disease factor patterns within the 95th percentile mask after weighted factor compositions for each patient. The association between clinical improvements and brain features was constructed using the PLSR approach in the training set. The established model was used for the testing set to obtain the predicted clinical improvement in this individual. The above-mentioned process was repeated n times (n = 122), and then the predicted clinical improvement score was obtained for each participant. The permutation test was then used to determine the correlation between observed and predicted clinical improvements. The observed depressive symptom improvements (pretreatment-posttreatment)/pretreatment × 100%87 across longitudinal patients, were randomly relabeled 1000 times to test the null hypothesis that the observed improvements were not related to the predicted improvements other than by chance. Distribution of the test statistic was then obtained, and the Pperm was obtained using occupied null models (> 95th centile).

In addition, longitudinal MDD patients were partitioned into three scales based on treatment outcomes, to determine relationships between clinical outcomes and abnormal factors. The responder group was defined as clinical improvement [i.e., (pretreatment-posttreatment)/pretreatment × 100% depressive symptom score] higher than or equal to a 50% decrease from baseline depression scales to a trial endpoint. The partial responder group was defined as clinical improvement lower than 50% and higher than or equal to a 25% decrease from baseline depression scale scores. The non-responder group was defined as having clinical improvement lower than 25%. Longitudinal MDD patients were then assigned to three factors based on the “winner-take-all” rule of factor compositions.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Related Articles

Latent circuit inference from heterogeneous neural responses during cognitive tasks

Higher cortical areas carry a wide range of sensory, cognitive and motor signals mixed in heterogeneous responses of single neurons tuned to multiple task variables. Dimensionality reduction methods that rely on correlations between neural activity and task variables leave unknown how heterogeneous responses arise from connectivity to drive behavior. We develop the latent circuit model, a dimensionality reduction approach in which task variables interact via low-dimensional recurrent connectivity to produce behavioral output. We apply the latent circuit inference to recurrent neural networks trained to perform a context-dependent decision-making task and find a suppression mechanism in which contextual representations inhibit irrelevant sensory responses. We validate this mechanism by confirming the behavioral effects of patterned connectivity perturbations predicted by the latent circuit model. We find similar suppression of irrelevant sensory responses in the prefrontal cortex of monkeys performing the same task. We show that incorporating causal interactions among task variables is critical for identifying behaviorally relevant computations from neural response data.

Functional brain network dynamics mediate the relationship between female reproductive aging and interpersonal adversity

Premature reproductive aging is linked to heightened stress sensitivity and psychological maladjustment across the life course. However, the brain dynamics underlying this relationship are poorly understood. Here, to address this issue, we analyzed multimodal data from female participants in the Adolescent Brain and Cognitive Development (longitudinal, N = 441; aged 9–12 years) and Human Connectome-Aging (cross-sectional, N = 130; aged 36–60 years) studies. Age-specific intrinsic functional brain network dynamics mediated the link between reproductive aging and perceptions of greater interpersonal adversity. The adolescent profile overlapped areas of greater glutamatergic and dopaminergic receptor density, and the middle-aged profile was concentrated in visual, attentional and default mode networks. The two profiles showed opposite relationships with patterns of functional neural network variability and cortical atrophy observed in psychosis versus major depressive disorder. Our findings underscore the divergent patterns of brain aging linked to reproductive maturation versus senescence, which may explain developmentally specific vulnerabilities to distinct disorders.

Genome-wide analysis identifies novel shared loci between depression and white matter microstructure

Depression, a complex and heritable psychiatric disorder, is associated with alterations in white matter microstructure, yet their shared genetic basis remains largely unclear. Utilizing the largest available genome-wide association study (GWAS) datasets for depression (N = 674,452) and white matter microstructure (N = 33,224), assessed through diffusion tensor imaging metrics such as fractional anisotropy (FA) and mean diffusivity (MD), we employed linkage disequilibrium score regression method to estimate global genetic correlations, local analysis of [co]variant association approach to pinpoint genomic regions with local genetic correlations, and conjunctional false discovery rate analysis to identify shared variants. Our findings revealed that depression showed significant local genetic correlations with FA in 37 genomic regions and with MD in 59 regions, while global genetic correlations were weak. Variant-level analysis identified 78 distinct loci jointly associated with depression (25 novel loci) and FA (35 novel loci), and 41 distinct loci associated with depression (17 novel loci) and MD (25 novel loci). Further analyses showed that these shared loci exhibited both concordant and discordant effect directions between depression and white matter traits, as well as distinct yet overlapping hemispheric patterns in their genetic architecture. Enrichment analysis of these shared loci implicated biological processes related to metabolism and regulation. This study provides evidence of a mixed-direction shared genetic architecture between depression and white matter microstructure. The identification of specific loci and pathways offers potential insights for developing targeted interventions to improve white matter integrity and alleviate depressive symptoms.

Depression symptom-specific genetic associations in clinically diagnosed and proxy case Alzheimer’s disease

Depression is a risk factor for the later development of Alzheimer’s disease (AD), but evidence for the genetic relationship is mixed. Assessing depression symptom-specific genetic associations may better clarify this relationship. To address this, we conducted genome-wide meta-analysis (a genome-wide association study, GWAS) of the nine depression symptom items, plus their sum score, on the Patient Health Questionnaire (PHQ-9) (GWAS-equivalent N: 224,535–308,421) using data from UK Biobank, the GLAD study and PROTECT, identifying 37 genomic risk loci. Using six AD GWASs with varying proportions of clinical and proxy (family history) case ascertainment, we identified 20 significant genetic correlations with depression/depression symptoms. However, only one of these was identified with a clinical AD GWAS. Local genetic correlations were detected in 14 regions. No statistical colocalization was identified in these regions. However, the region of the transmembrane protein 106B gene (TMEM106B) showed colocalization between multiple depression phenotypes and both clinical-only and clinical + proxy AD. Mendelian randomization and polygenic risk score analyses did not yield significant results after multiple testing correction in either direction. Our findings do not demonstrate a causal role of depression/depression symptoms on AD and suggest that previous evidence of genetic overlap between depression and AD may be driven by the inclusion of family history-based proxy cases/controls. However, colocalization at TMEM106B warrants further investigation.

Person-centered analyses reveal that developmental adversity at moderate levels and neural threat/safety discrimination are associated with lower anxiety in early adulthood

Parsing heterogeneity in the nature of adversity exposure and neurobiological functioning may facilitate better understanding of how adversity shapes individual variation in risk for and resilience against anxiety. One putative mechanism linking adversity exposure with anxiety is disrupted threat and safety learning. Here, we applied a person-centered approach (latent profile analysis) to characterize patterns of adversity exposure at specific developmental stages and threat/safety discrimination in corticolimbic circuitry in 120 young adults. We then compared how the resultant profiles differed in anxiety symptoms. Three latent profiles emerged: (1) a group with lower lifetime adversity, higher neural activation to threat, and lower neural activation to safety; (2) a group with moderate adversity during middle childhood and adolescence, lower neural activation to threat, and higher neural activation to safety; and (3) a group with higher lifetime adversity exposure and minimal neural activation to both threat and safety. Individuals in the second profile had lower anxiety than the other profiles. These findings demonstrate how variability in within-person combinations of adversity exposure and neural threat/safety discrimination can differentially relate to anxiety, and suggest that for some individuals, moderate adversity exposure during middle childhood and adolescence could be associated with processes that foster resilience to future anxiety.

Responses

Your email address will not be published. Required fields are marked *