Rehabilitative compensatory mechanism of hierarchical subnetworks in major depressive disorder: A longitudinal study across multi-sites

Abstract

Background:

Brain structural connectome comprise of a minority of efficiently interconnected rich club nodes that are regarded as ‘high-order regions’. The remission of major depressive disorder (MDD) in response to selective serotonin reuptake inhibitor (SSRI) treatment could be investigated by the hierarchical structural connectomes’ alterations of subnetworks.

Methods:

Fifty-five MDD patients who achieved remission underwent diffusion tensors imaging (DTI) scanning from 3 cohorts before and after 8-weeks antidepressant treatment. Five hierarchical subnetworks namely, rich, local, feeder, rich-feeder and feeder-local, were constructed according to the different combinations of connections and nodes as defined by rich club architecture. The critical treatment-related subnetwork pattern was explored by multivariate pattern analysis with support vector machine to differ the pre-/post-treatment patients. Then, relationships between graph metrics of discriminative subnetworks/ nodes and clinical variables were further explored.

Results:

The feeder-local subnetwork presented the most discriminative power in differing pre-/post- treatment patients, while the rich-feeder subnetwork had the highest discriminative power when comparing pre-treatment patients and controls. Furthermore, based on the feeder connection, which indicates the information transmission between the core and non-core architectures of brain networks, its topological measures were found to be significantly correlated with the reduction rate of 17-item Hamilton Rating Scale for Depression.

Conclusion:

Although pathological lesion on MDD relied on abnormal core organization, disease remission was association with the compensation from non-core organization. These results suggested that the dysfunctions arising from hierarchical subnetworks are compensated by increased information interactions between core brain regions and functionally diverse regions.

1. Introduction

Major depressive disorder (MDD) affects over 300 million individuals worldwide [1] and is associated with high rates of relapse [2]. While first episode MDD patients have a 40–60% chance of recurrence, second episode sufferers have an approximate 60% of future relapse and individuals above three episodes have a recurrence rate as high as 90% [3]. Selective serotonin reuptake inhibitors (SSRI) are first-line antidepressants which are effective in the treatment of both first episode and recurrent MDD patients [4]. Although SSRI generally have fewer side effects and work by increasing serotonin levels in the brain, their effects on neuroanatomical and functional changes have yet been fully understood [5].

Diffusion tensors imaging (DTI) is a promising neuroimaging technique that has the ability to map microstructural changes in white matter tracts and characterize the diffusion of water as a function of spatial location in response to antidepressant medications [6]. A number of DTI studies have examined the effects of antidepressants on abnormal brain regions or neural circuits and reported some region-specific treatment-induced alterations [7–9]. Following antidepressant pharmacotherapy, fractional anisotropy (FA) values indicated that the connectivity of the anterior cingulate and stria terminalis tracts were altered, suggesting that the anterior cingulate–limbic white matter tracts is a useful predictor of antidepressant treatment outcome in MDD [7]. In addition, remitted MDD patients presented with a volume increase in the left amygdala compared to healthy controls [10]. Treatment-resistant depression patients were compared to MDD patients before treatment and significant reductions in the FA of the bilateral hippocampus were observed [9].

On a network level, successful treatment with escitalopram or cognitive behavioral therapy was linked to the modulation of resting-state brain activity or metabolic changes in brain areas within the fronto-limbic circuit in MDD patients [5, 11]. There was a frontoparietal network response inhibition in MDD patients during remission with SSRI [12]. It was suggested that antidepressant medication was associated with the modular reconfiguration of hubs within the fronto-limbic circuit [13]. Furthermore, structural connectomic alterations were found between nodes of the default mode network (DMN) and fronto-thalamo-caudate brain regions in depressed patients [8]. Within the subnetworks of the DMN, a selective effect has been reported after the administration of antidepressants. Particularly, abnormal functional connectivity persisted within the anterior DMN subnetwork whereas the posterior DMN subnetwork was normalized after antidepressant treatment in depressed individuals [14]. The finding insinuated that the different subnetworks could perform different roles in the remission of MDD. However, the brain hierarchical modules, information communication and integration between modules are still poorly understood upon the administration of SSRI. The remission of MDD after the taking antidepressants could provide valuable information on numerous aspects of hierarchical subnetworks.

Rich club architecture, which is an organization of group of interconnecting nodes, plays an essential role in ‘high-order’ topology of brain networks and hierarchical subnetworks [15, 16]. The connections of rich club could be analyzed from densely interconnected ‘hub complex’ or ‘top hierarchical module’ [17]. The partitions that accounts for hierarchical subnetworks could be constructed by integrating different types of nodes and edges. In this way, the organization of rich clubs could be used to define hierarchical subnetworks and subsequently, the integration of information could be presented between distinct levels of brain networks. So far, most studies that are related to rich club organization focused on the disturbance or abnormality during disease stages and the alterations concerning rich club organizations during the remission period were overlooked. It has been ascertained time and again that the rich club networks were altered in various psychiatric diseases including schizophrenia [18, 19], and Alzheimer’s disease [20, 21], bipolar disorder [22], late-life depression [23]. However, there are few MDD studies which explored rich club organization after SSRI treatment.

In the current study, we utilized DTI to evaluate the effect of antidepressant treatments on hierarchical subnetwork configurations based on the rich club architecture. Our main objective was to determine how brain hierarchical subnetworks change after antidepressant treatment and examine whether MDD patients have abnormal rich club organization which is analogous to other psychiatric diseases. Therefore, five hierarchical subnetworks which focused on different levels of topological centrality and information integration were constructed by combining types of nodes and connections. The roles of these subnetworks were explored via the multivariate pattern analysis (MVPA) in the form of support vector machine. Afterwards, the patterns of pathological lesions and rehabilitation were differentiated so as to observe the hierarchical subnetworks along with their reconfiguration during SSRI pharmacotherapy in MDD patients.

2. Methods

2.1. Participants

The overall samples comprised of three datasets of participants, totally containing 165 DTI data from 55 MDD patients who underwent twice scans and 55 matched healthy controls from local community. 110 subjects, including 55 MDD patients were recruited from the Affiliated Brain Hospital of Nanjing Medical University and Nanjing Drum Tower Hospital and 55 matched healthy controls were enrolled through advertisements in the local community. Participants were classified into 3 cohorts namely, cohort A and B from the Affiliated Hospital of Nanjing Brain Hospital and cohort C from Nanjing Drum Tower Hospital.

Patients were diagnosed with MDD according to the criteria of the Diagnostic and Statistical Manual of Mental Disorders, 4th edition (DSM-IV) [24] and the severity of depression was assessed by the 17-item Hamilton Rating Scale for Depression (HAMD) [25]. Patients were included in this study according to the following criteria: (1) a total score of > 17 on the 17-item HAMD; (2) no comorbidity with other DSM-IV Axis-1 psychiatric disorders such as schizophrenia, substance abuse, obsessive compulsive disorders, and generalized anxiety; (3) no physical therapy for the last six months; (4) no use of psychotropic medication such as antidepressants, antipsychotics and benzodiazepines for the past two weeks. Subjects with the any of the following exclusion criteria were excepted from the study: (1) a history of alcohol or drug dependence; (2) serious medical illnesses such as organic brain disorders and severe somatic disease; (3) current pregnancy or breastfeeding; (4) inability to undergo magnetic resonance imaging (MRI) scanning. Healthy controls were screened using the Structured Clinical Interview for DSM‐IV Axis I disorders, research version, non‐patient edition, (SCID‐I/NP) and were omitted based on similar exclusion criteria. This research was undertaken in accordance with the ethical principle of the World Medical Association Declaration of Helsinki [29] and was authorized by the Research Ethics Review Board of each site. A complete description of the study was given and written informed consents were obtained from all subjects.

Patients from the three cohorts received an 8-week single-drug pharmacological treatment with SSRI. After the treatment, patients who achieved remission, defined as an HAMD scores < 7 [26] or > 50% reduction in HAMD scores [27], were enrolled for further scanning.

Cohort A consisted of 25 MDD patients who underwent 1.5 T MRI scanning. Of the initial 25 patients, 4 did not complete a second MRI scan, 1 did not achieve remission and 2 were excluded due to poor MRI signal quality. Therefore, a total of 18 participants were included in cohort A for analysis. Cohort B comprised 28 MDD patients who were submitted to 3 T MRI scanning. However, 2 patients refused to participate in second scans, 3 patients received a physical therapy due to their illness condition and 3 patients had excessive head movements during scanning. The remaining 20 patients in cohort B were included for further analysis. Cohort C included 26 drug-naive MDD patients were subjected to 3 T MRI scanning. Among them, 3 patients refused to undertake second scanning, 4 patients switched medications due to a poor response with SSRI and 2 patients had excessive head movements during scanning. Subsequently, a total of 17 patients remained in cohort C for analysis. Ultimately, 55 MDD patients and 55 healthy controls from the 3 cohorts were retained. Table 1 summarizes the demographic characteristics of subjects from the cohorts A, B and C.

Table 1 Demographic and clinical characteristics of participants.

Data are presented as the range of mean ± SD.

a The p value was obtained by two-sample two-tailed t test.

b The p value was obtained by two-tailed Pearson chi-square tests.

2.2. Treatment

The 55 patients who were finally included in the analysis received a single-drug treatment with SSRI. No systematic psychological therapies or stimulated intervention were performed. All patients were treated with escitalopram or sertraline according to a flexible dose antidepressant treatment protocol based on patient’s condition. The dose of escitalopram was gradually increased to 10–20 mg/day within a week and continued at this dose until finishing the study. And sertraline was started at 100 mg/day and could be increased to maximum of 200 mg/day. Of 55 patients, 48 patients were on escitalopram (mean dose = 18.1 mg/day) and 7 patients were on sertraline (mean dose = 142.7 mg/day). Detailed dosage information in each cohort was available in the Table S1. Since there were a restricted number of patients who received sertraline, the statistic was not performed on sertraline dosage. There were no significant differences of escitalopram dose (F = 0.499, p = 0.61) and medication type (χ² = 3.683, p = 0.159) across three cohorts.

2.3. Data acquisition

MRI data from Cohort A were acquired with a GE Sigma 1.5 T MRI scanner. Parameters for T1-weight axial images were: repetition time/echo time (TR/TE) = 500/14 ms, flip angle = 15°, inversion time = 400 ms, matrix = 256 × 256, voxel size = 1 × 1×2 mm³ and field of view (FOV) = 240 mm × 240 mm. For DTI scanning, diffusion was measured along 25 non-collinear directions (b value = 1000s/mm²) and an additional image without diffusion weighting (i.e., b = 0 s/mm²) was determined with parameters including TR/TE = 10000/81.2 ms, FOV = 240 mm × 240 mm, matrix = 128 × 128, voxel size = 1.9 × 1.9 × 4mm³ and slice thickness = 4 mm without gap.

Imaging data from Cohort B were obtained using a Siemens Verio 3 T MRI scanner. Parameters for T1-weight axial images were set as follows: TR/TE = 1900/2.48 ms, flip angle = 9°, inversion time = 900 ms, matrix = 256 × 256, voxel size = 1 × 1×1 mm³ and FOV = 250 mm × 250 mm. Diffusion in DTI scanning was measured along 30 non-collinear directions (b value = 1000s/mm²), and an additional image without diffusion weighting (i.e., b = 0 s/ mm²) was computed along with TR/TE = 6600/93 ms, FOV = 240 mm × 240 mm, matrix = 128 × 128, slice thickness = 3 mm without gap, voxel size = 1.9 × 1.9 × 3 mm³ and flip angle = 90°.

MRI data from Cohort C were acquired via a Siemens Verio 3 T MRI scanner. Parameters for T1-weight axial images included TR/TE = 2300/3.01 ms, flip angle = 9°, inversion time = 900 ms, matrix = 256 × 256, voxel size = 1 × 1×1 mm³ and FOV = 256 mm × 240 mm. DTI scanning parameters were as follows: TR/TE = 6100 /93 ms, FOV = 240 mm × 240 mm, matrix = 128 × 128, slice thickness = 2.8 mm without gap, voxel size = 1.9 × 1.9 × 2.8mm³ and flip angle = 90°.

2.4. Data preprocessing

Imaging data preprocessing were performed with the diffusion tool box of Functional Magnetic Resonance Imaging of the Brain (FMRIB) software library (FSL, http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/). The main procedures of preprocessing comprised four main parts namely, eddy current, motion artifacts correction, diffusion tensor estimation and FA calculation. Eddy current distortions and motion artifacts were corrected by applying a rigid body transformation from each diffusion-weighted image to the b₀ image. The diffusion tensor matrix was calculated according to Stejskal and Tanner equation [28]. Three eigenvalues and eigenvectors were obtained by diagonalization of the tensor matrix and the FA map was computed. Each DTI image was registered to the T1- weighted image and then to MNI-152 space. The transformation matrix from the diffusion space to MNI space was calculated by the aforementioned processing steps and was stored for further use.

2.5. Methodology overview

In order to simplify the proposed framework, a schematic overview of all the procedures is presented in Fig. 1. A MVPA explored the treatment-related subnetworks. MVPA inputs were distinct subnetworks’ structural graph-related characters (See Fig. 1B). The subnetworks which yielded the strongest discriminative abilities in the comparison of each two groups of pre-treatment patients, post-treatment patients and healthy controls, were detected following MVPA (Fig. 1C).

Fig. 1. Schematic overview of MVPA based on hierarchical subnetworks. A: Structural networks constructed from DTI modality. B: Five hierarchical subnetworks delineated by combinations of rich/un-rich nodes and their connections. Rich club nodes and connections are shown in red, non-rich club nodes in grey, feeder connection in orange and local connection in yellow. C: SVM analysis with graph kernel to differentiate each two groups of the pre-treatment patients, post-treatment patients and healthy controls (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).

2.6. Hierarchical subnetwork configuration based on Rich Club

Brain network nodes were defined using Automatic Anatomical Label (AAL) template [29], in which the whole brain was parceled out into 90 regions. Each edge of the network was defined by the mean FA value of the connected streamline between two regions.

In previous studies, rich club was measured by normalizing its coefficient at an individual level [30] or selecting the top 10%–20% high-degree nodes on a group level [21, 31]. In the current study, rich club nodes were defined as the top 20% high-degree nodes on a group-averaged weighted matrix M_group over the entire healthy samples. The M_group was constructed based on the rationale that only connections which were presented in a percentage of the total samples, specified by group threshold, were considered as present. A value of 50% was selected as the optimal group threshold across the range of 30% to 90% with a step size of 10% (see Supplementary Section 2 and Figure S1). M_group connection was defined by selecting the connections that were present in at least 50% of the all controls from the three cohorts. Element values for M_group were calculated by averaging the FA values of the selected connections across all subjects while the remaining connections were set to zero.

According to rich club definition, the nodes were categorized into two types: rich club nodes and non-rich club nodes. Network connections were classified into three types: the rich club connections that link members of rich club nodes, the feeder connections that link the rich club nodes to non-rich club nodes, and the local connections that link the non-rich nodes to non-rich nodes. The rich club connections and local connections represent information communication among the nodes that belong to the same community, while the feeder connections denoted the information communication between two hierarchical levels of nodes. These nodes and connections were combined in different ways to construct five types of hierarchical subnetworks. They are (1) Rich network, constructed via the rich club nodes and rich club connections. (Fig. 1B (2)) (2) Rich-feeder subnetwork including the feeder connections and non-rich nodes together with the rich network. (Fig. 1B (1)) (3) Feeder subnetwork including the feeder connections and all nodes. (Fig. 1B (4)) (4) Feeder-local subnetwork consisting of the feeder connections, local connections and all nodes. (Fig. 1B (3)) (5) Local network consisting of the non-rich nodes and local connections (Fig. 1B (5)). Supplementary Figure S3 and Supplementary Table S2-3 showed the graphic representation of five subnetworks and details information.

2.7. Discriminative subnetwork and node/graphlet detection

Graphs of distinct subnetworks were used to distinguish pre-treatment patients, post-treatment patients and healthy controls. In the MVPA scheme, the graph kernel extracted the similarity of graphs’ isomorphism as features based on Weisfeiler-Lehman sequence of graphs. The calculation procedures were described in Section 3 of supplementary material.

Binary support vector machine(SVM) classifiers were utilized with graph kernel to differentiate pre-treatment patients, post-treatment patients and controls according to the type of subnetwork configuration. Leave-one-out cross-validation was performed, where one sample was left out as test data and the remaining was used to train the model. Furthermore, a nest 10-fold cross-validation was used to optimize the SVM parameters. The subnetwork that achieved the highest accuracy was designated as the discriminative subnetwork.

Discriminative nodes or graphlets (subgraphs consisting of a core nodes and its connected nodes) in each subnetwork for MVPA need to be further detected. According to the definition of graph kernels, each element of the feature mapping (see Supplementary Figure S2 and Section 3) represented information of the corresponding nodes or graphlets. The particular element of feature mapping was eliminated prior to shuffling subjects and repeating MVPA runs. Then, the labeling of two groups was randomly shuffled and MVPA over these shuffled subjects was applied 1000 runs. The null distribution of accuracy’s changes between complete feature set and feature set without this particular variable was obtained. According the distribution, the corresponding nodes or graphlets which have significant powers for group discrimination survived after false discovery rate (FDR) correction, were denoted as the discriminative nodes or graphlets.

2.8. Relationships between discriminative subnetworks/nodes and clinical variables

The relationship between discriminative measures and clinical variables were assessed at two distinct levels including network level and node level. Since we hypothesized that distinct subnetworks reflected information communication variability, global efficiency (GE) which was calculated as the reciprocal of average path length, was utilized to present information communication efficiency in the networks [32]. In order to evaluate whether discriminative nodes were associated with remission level, nodal degree was employed to provide information about graphs’ isomorphism and the communication capability was defined as the total number of connections [33]. Both GE and nodal degree were computed within each subnetwork which was presented as a weighted graph.

Since HAMD scores reduction rate reflected the change in depressive severity and the change in GE as well as degree was also determined for better comparison. The nodal degree change rate and GE change rate were defined as: represent nodal degree and GE of each patient before and after treatment.

As the relationship was implemented upon three cohorts, all measures thus were normalized using rescaling method which was so-called Min-Max scaling via the following equation: X’=(X-X _min)/(X _max-X _min). On the aspect of networks, HAMD scores were correlated with GE of each subnetwork in terms of absolute value and change rate, respectively. On the aspect of nodes, correlations between reduction rate of HAMD and of discriminative nodes were separately calculated 15 times (5 subnetworks × 3 cohorts). Multiple comparison for correlation analysis was employed using FDR correction. Due to the large number of correlation pairs in three cohorts, two evaluation indicators, the averaged correlation value together with a standard deviation (SD) and the percent of correlation pairs with corrected significance in each subnetwork, were taken into consideration, so that the subnetwork with a great relevance to remission level could be identified. Averaged correlation values together with SDs revealed overall distribution or dispersion of the correlation performance of subnetworks. The percent of correlation pairs with corrected significance was described as a ratio of the number of nodes with corrected significance and the number of all the discriminative nodes in the corresponding subnetwork.

All correlations were computed by Pearson’s correlation. Outliers were identified by bootstrapping the Mahalanobis distance of each observation from the bivariate mean. Three subjects with a Mahalanobis distance of > 6 were excluded from the study. Significant level was set at p < 0.05 for all correlation analyses and multiple comparison correction.

2.9. Reproducibility and validation

In order to validate the main results, rich club regions were amended and a different design of network characterizations were employed. A new frame of hierarchical subnetworks was created by changing the threshold of rich club selection and the discriminative performances on this new hierarchical system were reappraised. Additionally, four topological metrics based on graph metrics were calculated and were used as inputs in MVPA instead of the subnetwork structure. More details were presented in Section 4 of supplementary material.

3. Results

3.1. Rich club organization

Rich club regions were delineated according to the top 20% highest-degree (k > 15) regions from M_group, as listed in supplementary Table S2. In the current study, rich club organization comprised of 22 regions, including bilateral precuneus, putamen and thalamus which is in line with the previous rich club findings [15]. The Graphical representation of five subnetworks was presented in supplementary Figure S3 and regions in rich network and local network were presented in supplementary Table S2-3.

3.2. Discriminative subnetworks and their correlations with clinical variables

Fig. 2A shows three bar graphs depicting the discriminative performance between each two groups of pre-treatment patients, post-treatment patients and health controls across the three cohorts. When differentiating pre-treatment patients from post-treatment patients, the feeder-local subnetwork yielded the highest discriminative performance of 78% compared to four other types of hierarchical subnetworks. Overall, the rich-feeder subnetwork possessed the highest discriminative performance of 86% in pre-treatment patients vs healthy controls in contrast to the performance of other subnetworks configurations.

Interestingly, the three cohorts demonstrated consistent trends of subnetwork superiorities in MVPA, that is, the feeder-local subnetwork had the highest discriminative power in pre/post treatment group comparison, and rich-feeder subnetwork had the highest discriminative power in pre-treatment patients and healthy controls. Fig. 2 B represents a scatterplot for correlation analysis between subnetwork GE values or of feeder-local subnetwork was significantly correlated with reduction rate of HAMD scores(r = 0.278, p = 0.041).

3.3. Discriminative nodes/graphlets and their correlations with clinical variables

Distinct nodes/graphlets were detected through a significant power survived after FDR correction. In this study, no discriminative graphlets were found. This could be partly explained by the fact that the occurrence frequency of the compressed labels (see Supplementary Section 3) contributed little in the computation of two networks’ similarities as the compressed labels differed greatly between individuals.

Fig. 2. A: Bar charts of MVPA performance over the five types of hierarchical subnetworks among the three cohorts. B: Scatterplots illustrating significant association between subnetworks’ GE and clinical variables. The solid lines represent the best-fitting linear regression. Pre: Pre-treatment patients; Post: Post-treatment patients; HC: Healthy Controls.

The graphical representation of discriminative nodes in each subnetwork was presented in Fig. 3 . The Results of correlation between of discriminative nodes and reduction rate of HAMD were also depicted in Fig. 3. The bar chart represents the percent of correlation pairs with corrected significance and line chart summarizes the averaged correlation value together with SD. From Fig. 3, feeder subnetwork had the highest percent of correlation pairs with corrected significance. The overall distribution showed that the correlation performance in feeder subnetwork tend to be better than other four subnetworks, suggesting that the change in nodal degree in the feeder subnetwork would be more likely to be correlated with individual’s remission level. Moreover, in the feeder subnetwork, the corrected significance of the left dorsolateral prefrontal cortex (DLPFC) and right precuneus were consistently achieved across the three cohorts.

4. Discussion

In this study, hierarchical subnetworks were developed in order to explore the patterns of neural networks in MDD patients who achieved remission following SSRI treatment. We reported that neural networks exhibited a hierarchical organization in which each subnetwork played different roles in the disease’ lesion and recovery. In order to characterize distinct information transmission patterns with respect to the interaction between different levels of neural networks, five types of subnetwork configurations were constructed by the different combinations of rich club connections as well as feeder and local connections. The MVPA results on these networks suggested that the brain pathological mechanism and compensatory rehabilitation on structural brain networks are diverse. The feeder-local subnetwork was significant in pre-treatment patients compared with post-treatment patients, while the rich-feeder subnetwork played an important role in disease process.

Rich club regions, similar to hub nodes, are highly centralized cerebral regions which are usually extensively interconnected with other areas as a core architecture [34]. Our results suggested that the brain lesions associated with depression were characterized by disrupted rich club architectures. The brain regions in the rich-feeder subnetwork, comprising of the insula, dorsolateral superior frontal gyrus and precuneus, were strongly associated with the brain lesions present in MDD. These regions are mostly located in the default mode network (DMN), cognitive-emotional circuitry and fronto-parietal circuitry. In prior studies, DMN regions were central to prominent functional hubs [18, 35, 36] and were frequently associated with regional cerebral hyperperfusion [37] as well as cerebral glucose hypermetabolism [38, 39]. It was previously demonstrated that there was a change in DMN and its subnetworks along with their interactions in MDD patients [40]. In addition, extensively impaired information interactions of abnormal brain anatomical topological organizations in cognitive-emotional circuitry and fronto-parietal circuitry were present in MDD patients [41]. Therefore, the dysfunctions in the rich-feeder subnetwork could lead to an increased susceptibility to pathogenic processes in MDD [42].

Fig. 3. The graphical representation of discriminative nodes in each subnetwork and two evaluation indictors among five subnetworks. Correlations between reduction rate of HAMD and Δdegree of discriminative nodes were summarized among three cohorts. Rich club nodes are displayed in red and non-rich club nodes are shown in grey. The red line represents rich club connections, blue line depicts feeder connections and yellow line denotes local connections. The bar chart shows the percent of corrected correlation pairs used to compare the performance of the five subnetworks. The line chart indicates averaged correlation value with SD. Error bars represent standard deviations. Feeder subnetwork achieved the best performance according to both indictors. The Δdegree of DLPFC and precuneus in the feeder subnetwork were significantly correlated with remission levels across the three cohorts (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).

As an interesting expansion of the rich-club role in depression attack, our key finding is related to the important role of the feeder-local subnetwork on depression recovery in comparison with other four types of subnetworks. Our results indicated that the feeder-local subnetwork of MDD patients exhibited more significant treatment-related compensatory changes compared to other subnetworks. Furthermore, the global efficiency of feeder-local subnetwork was significantly correlated with MDD remission level. These findings provided a novel remission mechanism whereby the recovery of MDD prominently relies on the compensatory changes arising from information interaction between core and non-core architectures. A central network was created by rich club nodes which attracted, transmitted, altered and integrated global brain signals’ communications. The feeder connections used to receive and integrate information from functionally diverse (non-rich club) brain regions [16]. A possible explanation of the biological mechanism is that antidepressants induced specific neural modulation of pathologic circuits and allowed patients to recover from depressive symptoms [43]. Therefore, we suggested that recovery is organized to facilitate or promote the normalization of whole brain function by enhancing the information integration process of non-rich regions during treatment with antidepressants. The purpose is to obtain a high level of integration between functional subsystems and as a complementation to lesion regions.

In our findings, there were correlations with corrected significance between changes in HAMD scores and regional degree of the left DLPFC and right precuneus (rich nodes) in the feeder subnetwork. This suggested that connections from rich regions to peripheral regions i.e. feeder connections, played a crucial role in the rehabilitation of MDD patients. Prior structural and functional studies described the DLPFC and precuneus as biomarkers of treatment response [44–46]. The DLPFC, which is the executive component of brain connectome, is implicated in mood regulation along with cognitive function [47–50] and is substantially affected by antidepressant treatment. It was previously demonstrated that a lower FA in DLPFC was associated with poor antidepressant response in geriatric depression [51]. Moreover, there was an improvement in the structural integrity between the DLPFC and deeper limbic structures after treatment with repetitive transcranial magnetic stimulation [52], suggesting an inhibitory influence of DLPFC over limbic structure during affective control [53]. Previous functional neuroimaging studies revealed hyperconnectivities between the left DLPFC and posterior DMN in post-ECT MDD remitters [54], left DLPFC and amygdala in MDD patients who underwent psychotherapy [55] as well as increased connectivities between DLPFC and hypothalamus after eight weeks of SSRI treatment [56]. The changes from rich club node of DLPFC to other functional regions could be interpreted as a potential association of cognitive control over affective regions, suggesting that antidepressant treatment strengthened the top-down control ability over limbic regions and consequently contributed to the improvement of depressive symptoms [55, 57]. The precuneus, which is considered as a connecter hub, forms an initial rich club organization along with other subcortical regions and plays a critical role in between-module connectivity [15, 34, 58]. Since the topological reconfiguration of large-scale cerebral networks are affected by antidepressant medication [13, 59], between-module connections i.e. feeder connections from precuneus to other peripheral regions might be altered after antidepressant treatment. It was previously reported that there was an increased communicability between the right precuneus and non-rich regions of the left inferior parietal gyrus and left superior temporal gyrus in remitted MDD patients compared to pre-treatment patients [60]. Therefore, it could be inferred that there is a potential cognitive compensation of DPLFC and modulation of the precuneus, which acted as crucial between-module connecter in the remission of MDD with SSRI monotherapy.

In this study, subnetwork topological information and graph kernel were used as new similarity measures based on graph isomorphism in an attempt to investigate whether the graphs were topologically identical. While other common kernels, such as Gaussian radial basis function, mostly employed on vector-type data without considering the structural information regarding the graph itself, graph kernel preserved the topological information of distinct subnetwork. Given that different rich club selections and network characterization could question the reproducibility of these findings, we repeated the experiment with different threshold settings of rich club definitions and information encoding based on graph metric. The results were consistent across the three datasets, confirming the generalization of our findings.

5. Conclusion

Our results are consistent with previous studies which revealed that organization of hierarchical hub regions of the human brain network were vulnerable in depression. Moreover, after antidepressant treatment, we found that cerebral function could be compensated by improving information interactions between core and non-core architecture as well as communication within non-core architecture. It could be inferred that MDD remission with SSRI was achieved by enhancing the integration of information between functionally diverse brain regions and rich club regions. In conclusion, we validated that the ways of brain lesion and rehabilitation compensation in structural brain networks are diverse.

Conflict of interest

The authors state that they have no conflict of interest with the content of this article.

Funding

The work was supported by the National Natural Science Foundation of China [grant numbers: 81871066, 81571639]; Jiangsu Provincial Medical Innovation Team of the Project of Invigorating Health Care through Science, Technology and Education [grant number: CXTDC2016004]; Jiangsu Provincial key research and development program [grant number: BE2018609].