Gray matter atrophy is constrained by normal structural brain network architecture in depression

Abstract

Background

There is growing evidence that gray matter atrophy is constrained by normal brain network (or connectome) architecture in neuropsychiatric disorders. However, whether this finding holds true in individuals with depression remains unknown. In this study, we aimed to investigate the association between gray matter atrophy and normal connectome architecture at individual level in depression.

Methods

In this study, 297 patients with depression and 256 healthy controls (HCs) from two independent Chinese dataset were included: a discovery dataset (105 never-treated first-episode patients and matched 130 HCs) and a replication dataset (106 patients and matched 126 HCs). For each patient, individualized regional atrophy was assessed using normative model and brain regions whose structural connectome profiles in HCs most resembled the atrophy patterns were identified as putative epicenters using a backfoward stepwise regression analysis.

Results

In general, the structural connectome architecture of the identified disease epicenters significantly explained 44% (±16%) variance of gray matter atrophy. While patients with depression demonstrated tremendous interindividual variations in the number and distribution of disease epicenters, several disease epicenters with higher participation coefficient than randomly selected regions, including the hippocampus, thalamus, and medial frontal gyrus were significantly shared by depression. Other brain regions with strong structural connections to the disease epicenters exhibited greater vulnerability. In addition, the association between connectome and gray matter atrophy uncovered two distinct subgroups with different ages of onset.

Conclusions

These results suggest that gray matter atrophy is constrained by structural brain connectome and elucidate the possible pathological progression in depression.

Introduction

Modern neuroimaging studies reveal that patients with depression exhibit widespread gray matter volume (GMV) abnormalities almost spanning the whole brain (Depping et al., 2016; Schmaal et al., 2016, 2017). These findings raise an open question whether these structural abnormalities follow some pattern or are just irregularly distributed. Brain regions are not isolated but form a highly interconnected network, or ‘connectome’, that enables efficient communication among brain regions to support diverse behavioral and cognitive functions (Bullmore & Sporns, 2009). Recent studies propose that pathological propagation of brain disorders is also constrained by brain connectome architecture (Yau, Zeighami, Baker, Larcher, & Vainik, 2018; Zhou, Gennatas, Kramer, Miller, & Seeley, 2012).

An increasing number of studies demonstrate close association between pathological propagation and connectome topology in brain disorders. Neurodegenerative diseases, such as Parkinson's and Alzheimer's diseases, are hypothesized to involve prion-like spread of toxic misfolded protein aggregates through transneuronal transport (Frost & Diamond, 2010). In neurodegenerative diseases, through modern neuroimaging methods, a series of studies propose that pathological perturbations may begin with focal ‘epicenters’ whose normal connectome architecture derived from normal population resembles the pattern of brain tissue volume loss and then propagate to unaffected brain regions with strong connections to the putative disease epicenters (Brown et al., 2019; Yau et al., 2018; Zeighami et al., 2015; Zhou et al., 2012). These disease epicenters are preferentially anchored to brain ‘hub’ regions playing key roles in efficient communication among brain regions (Larivière et al., 2020). What is more, the normal connectome architecture can forecast longitudinal gray matter atrophy in neurodegenerative diseases (Brown et al., 2019). In addition to neurodegenerative diseases, recent works report that the association hold true in other conditions, such as normal aging and schizophrenia. For example, age-related thickness differences are found to be constrained by connectome architecture (Petersen et al., 2022). In schizophrenia, Wannan et al., demonstrate that cortical thickness reductions are more likely among brain regions with stronger structural covariance connections, providing a potential explanation for the irregularly distributed cortical thickness reductions (Wannan et al., 2019). Shafiei et al., further establish that structural tissue volume loss is constrained by structural and functional connectome topology with a putative disease epicenter in the anterior cingulate cortex (Shafiei et al., 2020). These studies elucidate how the pathological perturbations propagate from the putative disease epicenters to affected brain regions, deepening our understanding of widespread pathophysiological effects and facilitating targeted treatment to prevent pathological progression in these diseases. In depression, studies consistently report progressive structural abnormalities in cortical and subcortical brain regions (Han et al., 2021b; Morris et al., 2019; Schmaal et al., 2016, 2017; Verduijn et al., 2015). However, whether structural abnormalities are associated with connectome topology remains unknown in depression.

On the other hand, although neuroimaging studies have consistently identified widespread GMV abnormalities in depression, their findings are not uniform (Chen et al., 2016; Han et al., 2022d). The reason is that mental disorders including depression are heterogeneous disorders and patients with mental disorders exhibit high interindividual variation in structural abnormalities (Kessler et al., 2003; Yu et al., 2019). To handle with the neuroanatomical heterogeneity, a series of methods are proposed, such as revealing more homogeneous subgroups (Beijers, Wardenaar, van Loo, & Schoevers, 2019; Han, Xu, Guo, Fang, & Wei, 2022a, 2022b; Lynch, Gunning, & Liston, 2020) and even inferring individual-level structural differences (Han et al., 2022d; Marquand, Rezek, Buitelaar, & Beckmann, 2016). These studies suggest that the group-level neuroanatomical differences are not representative of most patients (Lv, Di Biase, Cash, & Cocchi, 2020; Sun et al., 2021). In this context, one the limitations of these studies investigating the association between brain connectome architecture and structural brain abnormalities is that most of them just focus on the association at group level (Shafiei et al., 2020; Wannan et al., 2019), whether their findings hold true at individual level remains unknown.

In this study, we investigated the association between healthy connectome architecture and gray matter atrophy at individual level in depression. Two independent datasets were included this study: a discovery dataset (105 never-treated first-episode patients with depression and matched 130 healthy controls [HCs]) and a replication dataset (106 patients with depression and matched 126 HCs). First, the gray matte volume (GMV) was measured with voxel based morphometry analysis (VBM) (Ashburner & Friston, 2000) and individualized differences were derived with normative model (Marquand et al., 2016) that could infer individual-level GMV differences by measuring its deviates from the healthy distribution (Marquand et al., 2016). Then, we identified the putative disease epicenter(s) for each patient and significantly shared disease epicenter(s) across patients in depression. According to the previous relevant studies, we hypothesized that these shared disease epicenter(s) were anchored to brain ‘hub’ regions. Third, we assumed that brain regions with stronger connections with the identified disease epicenter(s) would show greater vulnerability to disease.

Methods

Participants

Two independent datasets were included in this study: a discovery dataset and a replication dataset. All results reported in the main text were based on the discovery dataset and reproduced on the replication dataset. The discovery dataset included 105 never-treated first-episode patients with depression and matched 130 healthy controls (HCs) recruited from the first out-patient services of Department of Psychiatry, the First Affiliated Hospital of Zhengzhou University. The protocol to recruit the discovery dataset was approved by the research ethical committee of the First Affiliated Hospital of Zhengzhou University. The replication dataset included 106 patients with depression and matched 126 HCs. The protocol to recruit the replication dataset was approved by the research ethical committee of the University of Electronic Science and Technology of China. All patients were under depressive state. All study protocols were performed according to the Helsinki Declaration of 1975 and written informed consents were obtained from all participants before the experiment.

In the discovery dataset, patients were diagnosed according to Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV) for depression by one chief physician and one well-trained psychiatrist. The exclusion criteria included: comorbidity of other mental/psychotic disorders and a history of manic symptoms. The clinical states and severity of symptoms were evaluated using the 17-items Hamilton Depression scale (HAMD). Healthy controls (HCs) were recruited from the community through poster advertisement. None of them presented a history of serious medical or neuropsychiatric illness or a family history of major psychiatric or neurological illness in their first-degree relatives. Patients and HCs were all Han Chinese and right handedness. For all participants (both HCs and patients), the additional exclusion criteria included: (1) taking drugs, such as anesthesia, sleeping and analgesia in the past 1 month; (2) substance abuse; (3) a history of brain tumor, trauma, surgery, or other organic body disease; (4) suffering from cardiovascular diseases, diabetes, hypertension; (5) contraindications for MRI scanning and (6) other structural brain abnormalities revealed by MRI scan.

In the replication dataset, 106 patients with depression were recruited from the Clinical Hospital of Chengdu Brain Science Institute, University of Electronic Science and Technology of China. All patients were interviewed by two experienced psychiatrists using the Structured Clinical Interview for DSM-IV-TR-Patient Edition (SCID-P, 2/2001 revision) for depression. The 24-Item Hamilton Depression Scale was used to evaluate the clinical state. The exclusion criteria included: comorbidity with schizophrenia, mental retardation, personality disorder or anxiety, a history of loss of consciousness, substance abuse, and serious medical, neurological illness. All patients were treated with one of the selective serotonin and serotonin–norepinephrine reuptake inhibitors. HCs (n = 126) were recruited from the community through poster advertisements and interviewed using SCID (nonpatient edition). No healthy participant had a history of serious medical or neuropsychiatric illness and a family history of major psychiatric or neurological illness in their first-degree relatives.

Scan acquisition

Discovery dataset

The functional and T1-weighted anatomical images were acquired on 3-Tsela GE Discovery MR750 scanner (General Electric, Fairfield Connecticut, USA). The scanning parameters of T1-weighted anatomical images were as follow: repetition time = 8164 ms, inversion time = 900 ms, echo time = 3.18 ms, flip angle = 7 degrees, voxel size = 1 × 1 × 1 mm³, resolution matrix = 256 × 256, slices = 188, thickness = 1.0 mm. Functional images were obtained using an echo-planar imaging sequence. The parameters were as follow: TR/TE = 2000/40 ms, 32 slices, matrix size = 64 × 64, voxel size = 3.75 × 3.75 × 3.2 mm³, flip angle = 90°, slice thickness = 4.0 mm, gap = 0.5 mm, and total 180 volumes.

Replication dataset

The structural and functional images were acquired on a 3-Tesla GE Discovery MR750 scanner (General Electric, Fairfield Connecticut, USA). Structural T1-weighted images with 3D spoiled gradient echo scan sequence with the following parameters: TR/TE = 5.92/1.956 ms, voxel size = 1 mm × 1 mm × 1 mm, slice thickness = 1 mm, no gap, flip angle = 12°, matrix size = 256 × 256, and 156 slices. Functional images were obtained using an echo-planar imaging sequence with the following parameters: TR/TE = 2000/30 ms, 43 slices, matrix size = 64 × 64, voxel size = 3.75 × 3.75 × 3.2 mm³, flip angle = 90°, slice thickness = 3.2 mm, no gap, and total 255 volumes.

Data preprocessing

Functional images were preprocessed using SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/) and Data Processing Assistant for Resting-State fMRI package (DPABI, http://www.restfmri.net). Preprocessing steps included: removing the first 10 volumes, slice timing and realignment. Subjects were excluded if the translational and rotational displacement exceeded 3.0 mm or 3.0°. None of the participants were removed for this criterion. Then, functional images were then normalized to the standard EPI template (voxel size 3 × 3 × 3 mm³). Then images were smoothed using a 6 × 6 × 6 mm³ full width at half maximum (FWHM) Gaussian kernel and detrended to reduce low-frequency drift. Band-pass filtered (0.01–0.1 Hz) to remove high-frequency physiological noise. Nuisance covariates including Friston 24 motion parameters (Satterthwaite et al., 2012), white matter signal and cerebrospinal fluid signal were regressed out. Finally, the outliners were placed with the best estimate using a third-order spline fit to clean the portions of time course. Outliners were detected based on the median absolute deviation, as implemented in 3dDespike (http://afni.nimh.nih.gov/afni) (Allen et al., 2014). For each subject, the mean frame-wise displacement (FD) was calculated (Han et al., 2017, 2018). There was no significant difference in terms of mean FD between patients with depression and HCs with two sample t test (p values >0.05 for these two datasets).

Voxel based morphometry analysis

GMV was measured with VBM analysis (Ashburner & Friston, 2000). This procedure was done using the recommended pipelines with the Computational Anatomy Toolbox (CAT12, http://dbm.neuro.uni-jena.de/cat12/). We just briefly described the main steps here. The structural images were segmented into gray matter, white matter, and cerebrospinal fluid and normalized into Montreal Neurological Institute space (voxel size 1.5 mm × 1.5 mm × 1.5 mm). The structural images were nonlinearly modulated and smoothed with an 6 mm full width at half maximum Gaussian kernel (Ashburner, 2009; Han et al., 2021b) (Fig. 1a). The total intracranial volume (TIV) was also obtained for the following analysis. The Image Quality Rating (IQR) combing image noise contrast ratio, inhomogeneity contrast ratio and root mean square was recorded to assess images quality (Brown et al., 2019; Han et al., 2023a; Han et al., 2022a, b).

Figure 1. Study design schematic. (a) For each patient, gray matter volume (GMV) is measured using voxel based morphometry analysis. (b) and (c) The regional GMV differences (Z scores) at individual level are obtained using normative model. (d) The optimal model (structural v. functional, single v. multi) is determined. (e) The disease epicenters significantly shared by depression are identified. (f) The topological properties of the identified disease epicenters are examined. Disease exposure is calculated and correlated with regional Z scores to investigate the possible pathological progression from the disease epicenters to other unaffected brain regions.

Determination of individualized gray matter atrophy pattern

For each patient, we obtained personalized gray matter atrophy pattern for each patient with depression using normative model by evaluating its deviates from the normal distribution. Following the previous studies (Bethlehem, Seidlitz, Romero-Garcia, & Trakoshis, 2020; Shan et al., 2022), we trained a Gaussian process regression model to predict regional GMV from age and sex in HCs and then applied the trained model to infer regional GMV of patients. For each brain region, the Z-score was obtained by quantifying the deviation between the true GMV and the predicted one then normalized by the uncertainty of prediction (see the equation below).

$${\rm Zscore} = \displaystyle{{{\rm GM}{\rm V}_{{\rm predicted}}-{\rm GM}{\rm V}_{{\rm true}}} \over {\rm \sigma }}$$

where σ was the square root of variance (estimated from the Gaussian process regression). The positive Z-scores indicated gray matter atrophy in patients with depression compared with HCs and vice versa (Fig. 1b and 1c).

Healthy functional and structural covariance network generation

In previous studies, both structural network measured with structural covariance and functional network measured with functional connectivity were found to be associated with the tissue volume loss in schizophrenia (Shafiei et al., 2020; Wannan et al., 2019). Thus, we first investigated which one (structural network or functional network) better recapitulated gray matter atrophy in depression.

Functional connectivity networks were generated by computing pairwise Pearson's correlations between the average time series of 246 brain regions followed by Fisher'z transformation in healthy population (http://www.brainnetome.org/) (Dong et al., 2022; Fan et al., 2016). This brain atlas was chosen for its joint validation using both functional and structural anatomy and connectivity (Fan et al., 2016). Four brain regions were excluded for low SNR (average SNR <100). As a result, we obtained a 242 × 242 × N functional connectivity matrix for healthy group (N was the number of HCs).

Structural covariance network was obtained by correlating pairwise Pearson's correlations between the average GMV of 246 brain regions in healthy population (Han et al., 2023c). Age, sex and education level were first regressed from GMV (Wannan et al., 2019). The obtained correlation coefficients were transformed to Fisher'z scores to improve normality. To be consistent with functional connectivity, brain regions with low SNR identified before were also excluded. This yielded a structural covariance matrix (242 × 242) for healthy group.

Identification of disease epicenters for each patient

According to previous studies, disease epicenter was the brain region whose connectivity map in healthy group most resembled the regional gray matter atrophy map (Brown et al., 2019; Yau et al., 2018; Zeighami et al., 2015; Zhou et al., 2012). Specifically, the brain region whose connectivity map showed the highest correlation with the brain atrophy pattern (the epicenter goodness of fit score, GOF) was selected as epicenter (Brown et al., 2019). Although relevant researches typically hypothesized that there was one single epicenter (Brown et al., 2019; Shafiei et al., 2020; Zeighami et al., 2015), this hypothesis was untested in depression. We investigated if there was one or multiple disease epicenters (single-epicenter model or multiple-epicenter model) in depression (Fig. 1d).

The following procedures were performed based on structural and functional network respectively. For the single-epicenter model, the disease epicenter was defined as the brain region with the highest GOF. Specifically, for functional network, we performed one-sample t test to obtain a single group-level statistical parametric (t-statistic, 242 × 242) matrix. For structural connectivity, we used the 242 × 242 (Fisher'Z) structural covariance network. For each brain region and in each patient, its GOF was the Pearson's correlation between each row of connectivity map (t-statistic for functional network or Fisher'Z for structural network, 242 × 1) and the regional gray matter atrophy (Z scores, 242 × 1). For the multiple-epicenter model, we performed a backfoward stepwise regression analysis to explore the relationship between connectivity maps (structural or functional network) and regional gray matter atrophy pattern (Z scores, 242 × 1). Stepwise regression could automatically identify important predictive variables eliminating the variables which are marginal important and cause multicollinearity in consideration of the highly correlation between regional connectome patterns (Mayblyum et al., 2021). A set of brain regions were selected as putative disease epicenters. The performance of these models was evaluated using a linear regression where connectivity maps of the identified disease epicenter(s) were independent variables and gray matter atrophy pattern was dependent variable. The explained variance (R²) and its significance (p value) were recorded. The performance of these models was compared using paired t test.

For the identified disease epicenter(s), we investigated whether they were significantly shared by patients with depression above chance using permutation testing (1000 times). For each identified disease epicenter, we counted how many patients (N_true) shared it. For each permutation, a random network preserving the degree, strength weight and distributions of true one (structural or functional network) was constructed (Rubinov & Sporns, 2011) and disease epicenter(s) was identified based on the random network. Similarly, we counted the number of patients (N_rand) sharing it. The p value was defined as the ratio of N_rand larger than N_true. Finally, disease epicenter(s) whose p value (corrected for Bonferroni) was less than 0.05 was considered to be significantly shared by patients with depression (Fig. 1e).

Topological properties of disease epicenters

Then, we examined the topological properties of the identified disease epicenter(s). Previous studies had established that brain ‘hubs’ were vulnerable to disease and propagated pathology to unaffected brain regions (Fornito, Zalesky, & Breakspear, 2015; Zhou et al., 2012). To test this hypothesis in depression, we examined whether the identified disease epicenter(s) had higher degrees and participation coefficients than other brain regions (Fig. 1f).

The degree of a node in the weighted graph was defined as the sum of its connection strengths (Rubinov, Ypma, Watson, & Bullmore, 2015):

$$D_i = \mathop \sum \limits_{\,j\in N_i}^n e_{ij}$$

where j was one of the connected nodes of node i (N _i) and e _ij was the connection strength between node i and node j.

The participation coefficient (PC) of node i was defined as:

$$PC_i = 1-\mathop \sum \limits_{\,j\in N_i}^n \left({\displaystyle{{D_{im}} \over {D_i}}} \right)^2$$

Where m was one of the brain network M, D _im was the sum of all connection strengths between i and network m and D _i was the degree as defined above (Yeo et al., 2011). PC of zero indicated a node that only connects with other nodes in its own module, whereas nodes with connections that were uniformly distributed across all modules have a PC of one. In this study, brain regions were grouped into Yeo-17 networks (Yeo et al., 2011). This procedure was conducted using brain connectivity toolbox (https://www.nitrc.org/projects/bct/).

To investigated whether shared disease epicenters had significantly higher degree and participation coefficient than other brain regions. The mean degree and PC of disease epicenters shared by patients with depression (identified in the previous step) were calculated and compared with those of randomly selected other regions with the same number. This procedure was repeated 1000 times and yielded a distribution of degree and PC of other brain regions.

Disease exposure

Then, we investigated the possible pathological progression from the identified disease epicenters to other unaffected brain regions. According to previous studies (Brown et al., 2019; Yau et al., 2018; Zeighami et al., 2015; Zhou et al., 2012), we hypothesized that regions with strong connections to disease epicenters would exhibit greater vulnerability defined by atrophy severity (Fig. 1f). For each patient and each region, the disease exposure was calculated (Yau et al., 2018).

The disease exposure was defined as the following equation:

$$Disease\;exposure( i ) = \mathop \sum \limits_j^N C_{ij} \times Atrophy( j ) $$

Where, C _ij was the connection strength (structural/functional connectivity strength) between unaffected brain region i with the disease epicenter j. Atrophy(j) was the Z-score of disease epicenter j. N was the number of the disease epicenters. Then, we calculated the Pearson's correlation coefficient between the disease exposure and atrophy severity (Z-scores) for each patient.

It was possible that the correlation between local deformation and the atrophy of connected neighbors was introduced by spatially proximal nodes that trivially exhibited greater co-deformation and greater connectivity. To rule out this possibility, we calculated ‘spatial disease exposure’ measure which was derived identically to the disease exposure except that the Euclidean distance was used instead of structural/functional connectivity strength.

Clustering analysis

Finally, we investigated whether the association between brain connectome architecture and gray matter atrophy pattern could help to uncover homogeneous subgroups in depression. As we found that structural multi-epicenter model best explained gray matter atrophy (see results), this procedure was performed based on structural network. Specifically, for each patient, the correlation coefficients (242 × 1) between each row of structural network (242 × 242) and regional gray matter atrophy (242 × 1, Z scores) were calculated and treated as features. The k-means clustering analysis was performed to uncover potential subgroups where the optimal number of clusters was determined (between 2 and 10) according to silhouette values. The k-means clustering was repeated 100 times to avoid local minima resulted from random in initialization of centroid positions (Allen et al., 2014). The squared Euclidean distance was used as the distance metric. We investigated whether demographic characteristics (including illness duration, age of onset, and symptom severity) differed among uncovered subgroups.

External validation and sensitivity analysis

To assess the reproducibility of our results, the main procedures above were performed on an independent replication dataset. Notably, in the replication dataset, patients had been treated with antidepressants and had different age distribution. As factors, including IQR and educational level exhibited significant differences in one of datasets (see results), we also investigated whether they affected the performance of models (R²) and association between disease exposure and gray matter atrophy using Pearson's correlation analysis.

Results

Demographic characteristics

The clinical demographics were presented in Table 1. There was no significant difference in terms of age and sex between patients with depression and HCs in both datasets. In the discovery, patients had significantly lower IQR than HCs. In the replication dataset, patients had significantly lower educational level compared with their respective control group.

Table 1. Demographic and clinical characteristics

Note: HC, healthy control; HAMD, Hamilton rating scale for depression; ^aχ2 test; ^btwo sample t test; IQR, Image Quality Rating.

Gray matter atrophy is constrained by structural network of focal disease epicenters

The performance (the explained variance, R ², mean ± s.d.) of these four models was as follow: structural single-epicenter model (0.12 ± 0.07), structural multi-epicenter model (0.44 ± 0.16), functional single-epicenter model (0.22 ± 0.03) and functional multi-epicenter model (0.09 ± 0.05). The structural multiple-epicenter model could significantly explain gray matter atrophy pattern for each patient (p < 0.001, Bonferroni corrected) and performed significantly better than other models (Fig. 2). The following results were all based on structural connectome.

Figure 2. The performance of models. The performance of models (explained R ²) are compared using paired t test (all p values <0.001). Note, SS, structural single-epicenter model; SM, structural multi-epicenter model; FS, functional single-epicenter model; FM, functional multi-epicenter model; SM-maximal atrophy, structural multi-epicenter model using connectivity maps of most atrophied brain regions.

Another question was that whether the identified disease epicenters were also the most atrophied brain regions, like the case in neurodegenerative diseases (Zhou et al., 2012). To answer this question, we investigated whether the connectome architecture of the identified disease epicenters showed closer relationship with gray matter atrophy than that of brain regions with maximal atrophy. Specifically, for each patient, the top N (the number of the identified disease epicenters) atrophied brain regions (according to Z scores) were selected. A linear regression was built where connectivity maps of brain regions with maximal atrophy were independent variables and gray matter atrophy pattern was dependent variables and the explained variance (R ²) was also recorded. We found that the connectome architecture of the identified disease epicenters recapitulated the gray matter atrophy better than brain regions with maximal atrophy did (for R ², paired t = 33.65, p < 0.001, Fig. 2).

Using structural multiple-epicenter model, we identified several disease epicenters for each patient. On average, patients with depression had 8.45 (± 2.97) disease epicenters (Fig. 3a). These disease epicenters exhibited tremendous interindividual variations in terms of number and distribution of disease epicenters. However, some of them were significantly shared by patients, such as the thalamus, hippocampus, striatum, medial frontal gyrus and inferior frontal gyrus (permutation testing p < 0.05, Bonferroni corrected) (orange nodes in Fig. 3b). The details were included in online Supplementary Table S1.

Figure 3. Interindividual variation in the distribution and number of disease epicenters and topological characteristics of disease epicenters significantly shared by depression. (a) Interindividual variation in the distribution and number of disease epicenters. (b) The distribution of disease epicenters. The node size represents the number of patients sharing this disease epicenter. The orange nodes represent disease epicenters significantly shared by depression (permutation testing p < 0.05, Bonferroni corrected). (c) Disease epicenters shared by depression exhibit higher participation coefficient than randomly selected other brain regions (p < 0.001). The orange line represents the mean PC of shared disease epicenters and the histogram represents the distribution of mean PC of randomly selected other brain regions. (d) Interindividual variation in the association between gray matter atrophy and disease exposure.

After identifying shared disease epicenters, we next sought to their topological characteristics. We found that these shared disease epicenters had significantly higher PC than other brain regions (permutation testing p < 0.001, Fig. 3c). However, there was no significant difference between these shared disease epicenters and other brain regions in terms of degree (permutation testing p > 0.05). The corresponding results of the replication dataset were included in online Supplementary Figure S1.

Brain regions with higher disease exposure show greater vulnerability to disease

On average, the regional disease exposure significantly correlated with gray matter atrophy pattern in patients with depression (the Pearson's correlation R, 0.50 ± 0.12, all p values <0.05, FDR corrected, Fig. 3d). To further rule out the effect of spatial proximity, we also defined the ‘spatial disease exposure’ and calculated its correlation with atrophy pattern. There was no significant correlation between ‘spatial disease exposure’ values and regional gray matter atrophy in depression (uncorrected p > 0.05 for each patient). To ensure that the correlation was not driven by nodes with high degree, the Pearson's correlation between degree distribution and averaged regional atrophy pattern across patients was also calculated. As a result, the correlation was not significant (the Pearson's correlation R = 0.06, p = 0.366).

The association between brain connectome architecture and gray matter atrophy uncovers homogeneous subgroups

Finally, we sought to investigate whether the association between connectome architecture and gray matter atrophy pattern could help to uncover homogeneous depression subgroups. As a result, patients with depression were optimally clustered into two subgroups (the corresponding silhouette values were drawn in Fig. 4A). Subgroup 1 (N = 43) significantly shared disease epicenters, such as thalamus, superior frontal gyrus, anterior cingulate cortex and parahippocampal gyrus (online Supplementary Table S2) while subgroup 2 (N = 62) shared disease epicenters, such as thalamus, superior frontal gyrus, anterior cingulate cortex, striatum and hippocampus (online Supplementary Table S3). The distribution of the disease epicenter was included in Fig. 4B. These two subgroups exhibited significant differences in terms of age of onset (two sample t = 2.48, FDR corrected p = 0.044, cohen’d = 0.49, Fig. 4C). In the replication dataset, patients were also divided into two subgroups (online Supplementary Figure S2).

Figure 4. The association between brain connectome architecture and gray matter atrophy uncovers two homogeneous subgroups. (a) The silhouette values for number of subgroups. (b) The distribution of disease epicenters for each subgroup. The node size represents the number of patients sharing this disease epicenter. The orange nodes represent disease epicenters significantly shared by depression (permutation testing p < 0.05, Bonferroni corrected). (c) These two subgroups exhibited significant differences in terms of age of onset.

External validation and sensitive analysis results

Most of our results could be reproduced in an independent dataset where patients had patients had been treated with antidepressants and had different age distribution. As the number (N = 31) of patients with younger age of onset (< = 18) was limited, the differences in terms of age of onset between the identified 2 subgroups tended to be significant (uncorrected p = 0.081, online Supplementary Figure S2). There was no significant difference between the performance of the optimal model (R ²), or association between disease exposure and GMV, and factors, including IQR and education level (all uncorrected p values >0.05).

Discussion

This study provided new insights regarding how the connectome architecture shaped gray matter atrophy in depression. We had three main findings. First, the distributed gray matter atrophy could be significantly explained by structural covariance profiles of the putative disease epicenters with higher participation coefficient in depression. Second, brain regions with stronger structural connections with the identified disease epicenters showed greater vulnerability to depression. Third, the association between connectome architecture and gray matter atrophy uncovered two homogeneous subgroups with difference age of onset. What's more, these findings could be reproduced in patients treated with medicine in an independent replication dataset suggesting the robustness of our findings.

Patients with mental disorders including depression demonstrate distributed structural atrophy. Whether these affected brain regions are irregularly distributed and how the pathological perturbations progress from the initial targets to other affected brain regions are of great concern. In schizophrenia and neurodegenerative diseases, recent studies reveal that patterns of structural atrophy are highly structured and constrained by normal brain network architecture of putative disease epicenter (Brown et al., 2019; Shafiei et al., 2020; Wannan et al., 2019; Yau et al., 2018; Zeighami et al., 2015; Zhou et al., 2012). In this study, we investigated this association in depression for the first time. Unlike most of relevant studies where the authors just hypothesize that these is only one disease epicenter, we first explored whether the multi-epicenter model or single-epicenter model better recapitulated structural atrophy pattern in depression. We found that the structural atrophy pattern could be better explained by the multi-epicenter model than single-epicenter model, suggesting there are multiple potential disease epicenters in depression. Although relevant studies report that gray matter abnormalities are related to both structural covariance and functional network, which one better recapitulate structural atrophy pattern in depression remains unknown (Shafiei et al., 2020; Wannan et al., 2019). Our another finding was that gray matter atrophy was better explained by normal structural covariance network relative to functional network. Mental disorders are associated with disturbed neurodevelopment trajectories among brain regions (Han et al., 2023c; Han et al., 2022d; Kaiser, Andrews-Hanna, Wager, & Pizzagalli, 2015; Lima-Ojeda, Rupprecht, & Baghai, 2018; Moberget et al., 2018; Sotiras et al., 2017; Yun et al., 2020). Compared with widely used functional connectivity, structural covariance measures brain connectivity features on a larger time scale (Evans, 2013) and potentially describes the coordinated regional volumes between brain regions reflecting their common development/maturation trajectories (Alexander-Bloch, Giedd, & Bullmore, 2013; Han et al., 2023b; Wang et al., 2023; Yun, Jang, Kim, Jung, & Kwon, 2015). Further analysis revealed that the identified disease epicenters showed significantly higher participation coefficient than other brain regions, explaining what topological characters made these regions to be disease epicenters. In our previous study, most of these brain regions, such as the hippocampus, thalamus and medial frontal gyrus, are found to exhibit GMV differences since the beginning of the disease (Han et al., 2021b). Another question is that whether the disease epicenters are also the most atrophied brain regions, like the case in neurodegenerative diseases (Zhou et al., 2012). The normal connectome architecture of the identified disease epicenters better explained brain atrophy than that of the most atrophied brain regions arguing that the disease epicenters were unnecessarily the most atrophied brain regions in depression. Taken together, our results suggest that pathological perturbation may preferentially target brain ‘hub’ regions probably due to their topological centrality, high metabolic demands and greater exposure to a toxic agent at the beginning of the disease in depression (Fornito et al., 2015).

In the neurodegenerative diseases, the pathological progression follows the normal network architecture as the transneuronal transport of toxic agents was associated with atrophy (Shafiei et al., 2020; Warren et al., 2013; Yau et al., 2018). More importantly, the network architecture can forecast the longitudinal atrophy in neurodegenerative diseases (Brown et al., 2019; Yau et al., 2018). In this study, we found that brain regions with strong structural covariance connections to the disease epicenters exhibited greater vulnerability. This finding provides a clue how the pathological perturbations progress from the initial targets to other affected brain regions. The mechanisms underlying this finding might be that disrupted interactions impede trophic support or normal neuronal signaling finally lead to atrophy of brain regions connected with disease epicenters (Fornito et al., 2015; Shafiei et al., 2020; Wannan et al., 2019).

A growing body of research recognize that depression was a highly heterogeneous disorder and the heterogeneity hamper neuroimaging studies reaching uniform findings in depression (Bondar, Caye, Chekroud, & Kieling, 2020; Drysdale, Grosenick, & Downar, 2017; Krishnan & Nestler, 2008). In neuroimaging studies, researchers acknowledge that patients with mental disorders exhibit tremendous inter-individual variation in morphological differences (Voineskos, Jacobs, & Ameis, 2020; Wolfers et al., 2018) and that the individual-level and group-level differences can be distinct and even opposite (Han et al., 2022d; Liu et al., 2021). However, most of previous studies investigate the association between brain atrophy and brain network architecture at group level (Shafiei et al., 2020; Wannan et al., 2019). Whether their findings held true at individual level remain understudied. For the first time, we investigated the association between morphological atrophy pattern and connectome architecture at individual level in depression. We elaborated this association at the individual level and showed that patients with depression demonstrated high inter-individual variation in the number and the distribution of disease epicenters. Using this individualized association, we uncovered two homogeneous subgroups with different age of onset, capturing the one of the sources of clinical heterogeneity. The age of onset is consistently recognized as an important potential confounder, leading to conflicting findings in depression (Chen et al., 2016; Han et al., 2021a; Han et al., 2021b; Herzog et al., 2021; Mai et al., 2021). For example, adult-onset depression is characterized by decreased volume in widespread brain regions (Schmaal et al., 2016) while adolescent-onset depression demonstrate increased volume in distributed brain regions, possibly reflecting an increased environmental sensitivity or delayed dendritic pruning (Blank & Meyer, 2022; Nickson et al., 2016). Our results show that the association between morphological atrophy pattern and connectome architecture also differ between adult- and adolescent-onset depression. These results provide new insights into the neuroanatomical heterogeneity and facilitate precision medicine for depression.

Our study has several limitations. First, we failed to record enough clinical data to deeply investigated whether and how the association between normal connectome architecture and gray matter atrophy underlay the clinical manifestation in depression. Second, only patients without other mental disorders comorbidities were recruited in this study. Future studies could recruit more heterogeneous samples to investigate the effect of comorbidity on these results. Third, in the replication dataset, all patients were taking medicine at the scan and we failed to record detailed information on medicine. Even our results suggested that findings derived from never-treated first-episode patients could be validated in patients who had taken antidepressants, we could not deeply investigate whether and how antidepressants affected these results.