Major depressive disorder on a neuromorphic continuum

Introduction

Major depressive disorder (MDD) is a devastating worldwide mental disorder¹. Patients with MDD are clinically diagnosed according to the DSM-5-TR when a patient exhibits at least five designated polythetic symptoms². 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 mechanisms³. However, such a categorical model does not consider alternative latent effects⁴. Indeed, depression has been more recently conceptualized as an unbroken “spectrum” (also referred to as factors), rather than as a “category”⁵. 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 symptoms⁶, neuroimaging intermediate phenotypes^7,8, and genotypes⁹. 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 skills¹⁰. 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 base^{11,12,13,14,15,16}. Progress toward translating knowledge into clinical practice has remained challenging¹⁷. One possible explanation may involve the clustering analytical approaches used to separate categories¹⁸. Indeed, categorical models merely highlight common brain features, whereas they disregard continuous individual variations within subtypes¹⁹. 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 boundaries^20,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 considered^25,26,27. Both exo-phenotype and genotype findings have confirmed the continuous variations across individuals, leading to high variance, even within subtypes^28,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 patient³¹. LDA has been successfully used to decompose whole-brain structural or functional features into abnormal patterns in various brain disorders^{32,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 MDD³⁸, 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 volumes^39,40. MDD is associated with widespread brain structural changes⁴¹. Cortical thickness is thought to reflect the density and arrangement of neurons, neuroglia, and nerve fibers⁴². Furthermore, meta-analyses have reported cortical thickness and subcortical volumetric abnormalities in MDD^8,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 regions^7,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.

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)⁴⁵. 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 systems^46,47.

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

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

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 images⁴⁸. 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 R² = 0.88, F_{(19, 62)} = 32.59, P = 2.2e-43, Fig. 4a; Factor 2: Adjusted R² = 0.92, F_{(19, 62)} = 49.29, P = 4.8e-50, Fig. 4b; Factor 3: Adjusted R² = 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-HT_2A (P = 0.0002, FDR-correction). The norepinephrine (P = 0.00003, FDR-correction), 5-HTT (P = 0.0002, FDR-correction), and mGluR₅ (P = 0.0006, FDR-correction) were primarily associated with Factor 2 distribution. For Factor 3, the serotonin system, including 5-HT_2A (P = 0.0006, FDR-correction), 5-HT₄ (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.**

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

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, P_perm = 0.01; Fig. 6a, left panel; pharmacological treatments: R = 0.33, P_perm = 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: χ² = 0.06, P = 0.97), and in pharmacological treatments (age: F = 0.35, P = 0.71; sex: χ² = 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.**

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 abnormalities⁴¹. However, many symptom domains manifest as continuous variations across individuals²⁶, even within subtypes⁴⁹. 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 findings⁵⁰, 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 regions⁵¹. 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 efficacies⁵². Other studies reported two converse MDD neuroimaging patterns related to default mode, limbic, sensorimotor, occipital, and subcortical regions^12,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 cortex^55,56. Specifically, our multivariate analyses found that the Factor 1 pattern was associated with the GABA_A receptor distribution map, crucial for maintaining neural balance⁵⁷. 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 regions⁵⁸. Recent studies have also suggested that cannabinoid receptor 1 was a promising therapeutic target for anti-depression and anti-anxiety treatments⁵⁹. The key brain regions in Factor 1, including orbitofrontal and occipital cortices, may become novel target areas for treating depressed patients^52,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 regulation⁶¹. Specifically, the CON was structurally and functionally altered in MDD patients⁶². 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 patients⁶³. 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 model⁶⁴, 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 illnesses⁶⁵. Furthermore, norepinephrine transporters, binding in the locus coeruleus, reflect the dynamics of norepinephrine and are decreased in postmortem tissues from patients with MDD⁶⁶, 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 controls^67,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 controversial⁶⁹, 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 depression⁷⁰. 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 treatments³. 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 patterns¹¹. 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 stages³⁴.

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 studies^71,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 samples⁷⁵ and was reduced in motion-affected group⁷⁶. 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 PCNtoolkit⁷⁷, 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 scans⁷⁵) as covariates⁷⁸.

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 function⁷⁹:

$$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|Factor_K), K = 1,2, …, K_max]³². 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 LDA³¹. The model used in the present study has been used to characterize latent disease factors in Alzheimer’s disease³⁴, epilepsy³², and autism study³³. 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 information^32,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 used⁸⁰, 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 study³². 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. ⁴⁸ and Markello et al. ⁸¹. These 19 receptors/transporters combined into nine neurotransmitter systems as follows⁸²: dopamine (D1, D2, DAT), norepinephrine (NET), serotonin (5-HT1_A, 5-HT1_B, 5-HT2, 5-HT₄, 5-HT₆, 5-HTT), acetylcholine (α₄β₂, M₁, VAChT), glutamate (NMDA, mGluR₅), GABA (GABA_A), histamine (H₃), opioid (MOR), and cannabinoid (CB₁). 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 factors^83,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 algorithm^33,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 features⁸⁶. 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 95^th 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%⁸⁷ 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 P_perm was obtained using occupied null models (> 95^th 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.