Resting-state connectivity biomarkers define neurophysiological subtypes of depression

Abstract

Biomarkers have transformed modern medicine but remain largely elusive in psychiatry, partly because there is a weak correspondence between diagnostic labels and their neurobiological substrates. Like other neuropsychiatric disorders, depression is not a unitary disease, but rather a heterogeneous syndrome that encompasses varied, co-occurring symptoms and divergent responses to treatment. By using functional magnetic resonance imaging (fMRI) in a large multisite sample (n = 1,188), we show here that patients with depression can be subdivided into four neurophysiological subtypes ('biotypes') defined by distinct patterns of dysfunctional connectivity in limbic and frontostriatal networks. Clustering patients on this basis enabled the development of diagnostic classifiers (biomarkers) with high (82-93%) sensitivity and specificity for depression subtypes in multisite validation (n = 711) and out-of-sample replication (n = 477) data sets. These biotypes cannot be differentiated solely on the basis of clinical features, but they are associated with differing clinical-symptom profiles. They also predict responsiveness to transcranial magnetic stimulation therapy (n = 154). Our results define novel subtypes of depression that transcend current diagnostic boundaries and may be useful for identifying the individuals who are most likely to benefit from targeted neurostimulation therapies.

Figure 1

Canonical correlation analysis (CCA) and hierarchical clustering define four connectivity-based biotypes of depression. (a) Data analysis schematic and workflow. After preprocessing, BOLD signal time series were extracted from 258 spherical regions of interest (ROIs) distributed across the cortex and subcortical structures. The schematics (top) show lateral (left) and medial (right) views of right-hemisphere ROIs projected onto an inflated cortical surface and colored by functional network (lower left). Left-hemisphere ROIs (data not shown) were similar. For each subject, whole-brain functional-connectivity matrices were generated by calculating pairwise BOLD signal correlations between all ROIs, as in this example of correlated signals (r² = 0.88) for DLPFC (solid line) and PPC (dashed line) nodes of the FPTC network in a representative subject. (b) Whole-brain, 258 × 258 functional-connectivity matrix averaged across all healthy controls (n = 378 subjects). z = Fischer transformed correlation coefficient. (c,d) CCA was used to define a low-dimensional representation of depression-related connectivity features and identified an “anhedonia-related” component (canonical variate; c) and an “anxiety-related” component (d), represented by linear combinations of connectivity features that were correlated with linear combinations of symptoms. The scatterplots in c and d illustrate the correlation between low-dimensional connectivity scores and low-dimensional clinical scores for the anhedonia-related (r² = 0.91) and anxiety-related components (r² = 0.95), respectively (P < 0.00001, n = 220 patients with depression). To the left of each scatterplot, clinical score loadings (i.e., the Pearson correlation coefficients between specific symptoms and the anhedonia- or anxiety-related clinical score (canonical variate)) are depicted for those symptoms with the strongest loadings (HAMD item ^#, indicated by numbers in superscript; for all loadings on all symptoms, see Supplementary Fig. 2). Below each scatterplot, connectivity score loadings are summarized by depicting the neuroanatomical distribution of the 25 ROIs (top 10%) that were most highly correlated with each component (summed across all significantly correlated connectivity features for a given ROI), colored by network, as in a. Projections to the medial wall map are for both left- and right-hemisphere ROIs. (e) Hierarchical clustering analysis. The height of each linkage in the dendrogram represents the distance between the clusters joined by that link. For reference, the dashed line denotes 20 times the mean distance between pairs of subjects within a cluster. For analyses of additional cluster solutions and further discussion, see Supplementary Figure 3. (f) Scatterplot for four clusters of subjects along dimensions of anhedonia- and anxiety-related connectivity. Gray data points indicate subjects with ambiguous cluster identities (edge cases, cluster silhouette values < 0; n = 15, or 6.8% of all subjects). ACC, anterior cingulate cortex; amyg, amygdala; antPFC, anterior prefrontal cortex; a.u., arbitrary units; AV, auditory/visual networks; CBL, cerebellum; COTC, cingulo-opercular task-control network; D/VAN, dorsal/ventral attention network; DLPFC, dorsolateral prefrontal cortex; DMN, default-mode network; DMPFC, dorsomedial prefrontal cortex; FPTC, frontoparietal task-control network; GP, globus pallidus; LIMB, limbic; MR, memory retrieval network; NAcc, nucleus accumbens; OFC, orbitofrontal cortex; PPC, posterior parietal cortex; precun, precuneus; sgACC, subgenual anterior cingulate cortex; SS1, primary somatosensory cortex; SN, salience network; SSM, somatosensory/motor networks; subC, subcortical; thal, thalamus; vHC, ventral hippocampus; VLPFC, ventrolateral prefrontal cortex; VMPFC, ventromedial prefrontal cortex; vStr, ventral striatum; n.s., not significant. See Supplementary Table 4 for MNI coordinates for ROIs in b and c.

Figure 2

Connectivity biomarkers define depression biotypes with distinct clinical profiles. (a) Neuroanatomical distribution of the 25 ROIs (top 10%) with the most abnormal connectivity features shared by all four biotypes (summed across all connectivity features for a given ROI), identified using Wilcoxon rank–sum tests to test for connectivity features that were significantly abnormal in all four biotypes relative to healthy controls in data set 1 (n = 378). ROIs are colored by network, as in Figure 1a. (b) Heat maps depicting a pattern of abnormal connectivity (P < 0.05, false-discovery rate (FDR) corrected) shared by all four biotypes for the top 50 most abnormal ROIs, colored on the basis of Wilcoxon rank–sum tests comparing patients and controls, as in a. Warm colors represent increase and cool colors decrease in depression as compared to controls. (c) Correlations (r = 0.72–0.82, ***P < 0.001, Spearman) between shared abnormal connectivity features (as indexed by the first principal component (PC) of the features depicted in b and the severity of the core depressive symptoms. Insets depict the prevalence of each symptom. Symptom severity measures are z-scored with respect to controls and plotted as the mean for each quartile, ± s.e.m. (d) Neuroanatomical distribution of dysfunctional connectivity features that differed by biotype, as identified by Kruskal–Wallis analysis of variance (ANOVA) (P < 0.05, FDR corrected), summarized for the 50 ROIs (top ~20%) with the most biotype-specific connectivity features (i.e., the 50 ROIs with the largest test statistic summed across all connectivity features, showing cluster specificity at a threshold of P < 0.05, FDR corrected). Nodes (ROIs) are colored to indicate the biotype with the most abnormal connectivity features and scaled to indicate how many connectivity features exhibited significant effects of biotype. (e) Heat maps depicting biotype-specific patterns of abnormal connectivity for the functional nodes illustrated in d, plus selected limbic areas, colored as in b. Green boxes highlight corresponding areas in each matrix discussed in the main text. (f) Biotype-specific clinical profiles for the six depressive symptoms that varied most significantly by cluster (P < 0.005, Kruskal–Wallis ANOVA). Symptom severities (HAMD) are z-scored with respect to the mean for all patients in the cluster-discovery set. See Supplementary Figure 4 for all 17 HAMD items and for replication in data from subjects left out of the cluster-discovery set. (g) Boxplot of biotype differences in overall depression severity (total HAMD score), in which boxes denote the median and interquartile range (IQR) and whiskers the minimum and maximum values. In f and g, asterisk (*) indicates significant difference from mean symptom severity rating for all patients (z = 0) at P < 0.05; error bars depict s.e.m.; n.s., not significant. Aud, auditory cortex; HC, hippocampus; lat PFC, lateral prefrontal cortex; lat OFC, lateral orbitofrontal cortex; MTG, middle temporal gyrus; PHC, parahippocampal cortex; PCC, posterior cingulate cortex; SSM, primary sensorimotor cortex (M1 or S1); STG, superior temporal gyrus; vis, visual cortex. Other abbreviations are as in Figure 1. See Supplementary Table 5 for Montreal Neurological Institute coordinates for ROIs in a and d.

Figure 3

Functional connectivity biomarkers for diagnosing neurophysiological biotypes of depression. (a) Data analysis schematic and workflow (Online Methods for additional details). (b) Optimization of diagnostic-classifier performance (accuracy) across the indicated combinations of methods for parcellation, clustering and classification. *P < 0.005, as estimated by permutation testing (Online Methods). Double asterisk (**) indicate the best performing protocol for parcellation, clustering and classification, and the focus of all subsequent analyses. (c–f) The neuroanatomical locations of the nodes with the most discriminating connectivity features are illustrated for each biotype for the four-cluster solution denoted by the double asterisk in b, colored and scaled by summing the results of Wilcoxon rank–sum tests of patients as compared to controls across all connectivity features associated with that node. Red represents increased and blue decreased functional connectivity in depression. (g) Sensitivity and specificity by biotype for the most successful classifiers identified in b (**). Error bars depict 95% confidence interval for the mean accuracy across all iterations of leave-one-out cross-validation. (h) Reproducibility of cluster assignments in a second fMRI scan (n = 50) obtained 4–5 weeks after the initial scan (χ² = 112.7, P < 0.00001). (i) Classifier performance in an independent, out-of-sample replication data set (n = 125 patients, 352 healthy controls). Cross-hatched bars depict classifier accuracy with more stringent data quality controls (Online Methods) and excluding equivocal classification outcomes (the 10% of subjects with the lowest absolute SVM classification scores). Error bars depict 95% confidence intervals.

Figure 4

Connectivity biomarkers predict differential antidepressant response to rTMS. (a) Differing response rates to repetitive transcranial magnetic stimulation (rTMS) of the dorsomedial prefrontal cortex across patient biotypes (clusters) in n = 124 subjects. Response rate indicates percentage of subjects showing at least a partial clinical response to rTMS (χ² = 25.7, P = 1.1 × 10⁻⁵), defined conventionally as >25% reduction in symptom severity by HAMD. Full response rates (>50% reduction by HAMD, cross-hatched bars) also varied by biotype (χ² = 22.9, P = 4.3 × 10⁻⁵). (b) Boxplot of percent improvement in depression severity by biotype (P = 1.79 × 10⁻⁶, Kruskal–Wallis ANOVA), in which boxes denote the median and interquartile range and whiskers the minimum and maximum up to 1.5 × the IQR, beyond which outliers are plotted individually. Percent improvement = total HAMD score before treatment – total HAMD score after treatment/total HAMD score before treatment. **P = 0.00001–0.002 (Mann–Whitney), indicating significantly increased versus biotypes 2–4; *P = 0.007 (Mann–Whitney), indicating significantly increased versus biotype 4. (c) Functional connectivity differences in the DMPFC stimulation target in treatment responders versus nonresponders (Wilcoxon rank–sum tests, thresholded at P < 0.005). Warm colors represent increased and cool colors decreased functional connectivity in treatment responders as compared to nonresponders. The 12 ROIs depicted here were located within 3 cm of the putative DMPFC target site, estimated in a previously published report to be located at Talairach coordinates, x = 0, y = +30, z = +30 (ref.). (d) The neuroanatomical distribution of the most discriminating connectivity features for the comparison of rTMS responders versus non-responders, summarized by illustrating the locations of the 25 (top 10%) most discriminating ROIs indexed by summing across all significantly discriminating connectivity features and colored by functional network as in Figure 1a. The red arrows denote the rTMS target site in the two (lower) medial panels. (e) Heat maps depicting differences in functional connectivity in patients who subsequently improved after receiving rTMS (n = 70), as compared to those who did not (n = 54). (f–i) Confusion matrices depicting the performance of classifiers trained to identify subsequent treatment responders on the basis of the most discriminating connectivity features (f), connectivity features plus biotype diagnosis (g), clinical symptoms alone (h) or connectivity features plus biotype diagnosis in an independent replication set (i, n = 30 patients with depression). NR, nonresponder; R, responder. (j) Summary of performance (overall accuracy) for classifiers in f–i. **significantly greater than clinical features alone (P < 0.001) and connectivity features alone (P = 0.003) by permutation testing; *P = 0.04 (significantly greater than clinical features alone by permutation testing). Cross-hatched bars depict classifier accuracy with more stringent data quality controls (Online Methods) and excluding equivocal classification outcomes (the 10% of subjects with the lowest absolute SVM classification scores). Error bars depict s.e.m. in a and 95% confidence intervals in j. All abbreviations as in Figures 1 and 2. See Supplementary Table 7 for MNI coordinates for ROIs in d.

Figure 5

Connectivity biomarkers of depression biotypes transcend diagnostic boundaries. (a) Abnormal connectivity features in patients with generalized anxiety disorder (GAD, n = 39) relative to healthy controls (n = 378). In this matrix depicting the 50 neuroanatomical nodes with the most significantly different connectivity features (Wilcoxon rank–sum tests, summed across all 258 features), elements in warm and cool colors depict connectivity features that are significantly increased or decreased in GAD, respectively. (b) 30.2% of connectivity features that were significantly abnormal in GAD (threshold of P < 0.001 versus controls, Wilcoxon) were also abnormal in depression (χ² = 5,457, P < 0.0001). (c) The neuroanatomical distribution of the most discriminating connectivity features for the comparison of GAD patients versus controls. The nodes are colored and scaled by summing across all significantly abnormal connectivity features associated with that node. Red represents increased and blue decreased functional connectivity in GAD. (d) Distribution of biotype diagnoses in patients with GAD. (e) No significant biotypes differences in anxiety symptom severity (P = 0.692; Kruskal–Wallis ANOVA). BAI, Beck anxiety inventory. (f,g) Significantly (P < 0.005, Kruskal–Wallis) elevated total depressive-symptom severity (f; BDI, Beck depression inventory) and anhedonia severity (g; BDI item 12) in GAD patients who tested positive for a depression biotype as compared to those who did not. *P < 0.01, ^†P = 0.064 in post hoc Mann–Whitney tests relative to “not depressed” group. (h) Distribution of biotype diagnoses in patients with schizophrenia (n = 41). Error bars depict s.e.m. throughout. All abbreviations as in Figures 1 and 2.