Resting-State Functional Connectivity Impairment in Patients with Major Depressive Episode

Aim: This study aims to develop new approaches to characterize brain networks to potentially contribute to a better understanding of mechanisms involved in depression. Method and subjects: We recruited 90 subjects: 49 healthy controls (HC) and 41 patients with a major depressive episode (MDE). All subjects underwent clinical evaluation and functional resting-state MRI. The data were processed investigating functional connectivity network measures across the two groups using Brain Connectivity Toolbox. The statistical inferences were developed at a functional network level, using a false discovery rate method. Linear discriminant analysis was used to differentiate between the two groups. Results and discussion: Significant differences in functional connectivity (FC) between depressed patients vs. healthy controls was demonstrated, with brain regions including the lingual gyrus, cerebellum, midcingulate cortex and thalamus more prominent in healthy subjects as compared to depression where the orbitofrontal cortex emerged as a key node. Linear discriminant analysis demonstrated that full-connectivity matrices were the most precise in differentiating between depression vs. health subjects. Conclusion: The study provides supportive evidence for impaired functional connectivity networks in MDE patients.


Introduction
Major depressive disorder (MDD) is a leading cause of disability worldwide [1]. However, its exact pathophysiological mechanism remains unclear. Symptoms of depressive episode include sadness, anhedonia, insomnia, restlessness, and suicidal thoughts [2]. Moreover, this condition is characterized by impairments in cognitive and emotional processing [3]. There is evidence suggesting that cognitive dysfunction could be seen not only in the acute phase but in remission as well, and it could even worsen the condition's outcome [4,5].
Brain networks can be defined as a set of regions that exhibit correlated activity in resting-state condition or during task performance [6,7]. In recent studies, researchers concentrated on building brain functional networks, as well as searching for abnormal communication between them in order to elucidate the pathophysiological mechanisms of MDD [8]. Different methods are used for constructing brain networks such as region of interest (ROI) analysis, seed-based analysis or independent component analysis (ICA) with the first two being preferred in hypothesis driven research; while the latter is an example of data driven approach [9]. However, for each of these methods there are strengths and weaknesses. In the ROI-based method the definition of the regions could affect functional connectivity patterns [10]. Moreover, seed-based analysis results depend much on the

Subjects
We recruited 90 subjects: 49 healthy controls (HC) and 41 patients with a major depressive episode (MDE) in the context of major depressive disorder (n = 35) or bipolar affective disorder (n = 6). Subjects from both groups were assessed by experienced psychiatrists using Mini International Neuropsychiatric Interview [27] and Montgomery-Åsberg Depression Rating Scale (MADRS) [28]. Subjects having a previous history of comorbid (for HC and patients, respectively) psychiatric conditions, autoimmune diseases, neurological diseases, history of head trauma, or any metal implants incompatible with the MRI were excluded. All participants provided a written consent form complying with the Declaration of Helsinki. The study was approved by the Medical University of Plovdiv Ethical Committee (2/19 April 2018).

Demographic and Clinical Characteristics
The two groups of subjects did not differ significantly in terms of mean age, sex, and level of education distributions. Expectedly, patients had significantly higher MADRS scores (see Table 1).

MR Scanning
The MR scanning procedure was performed on a 3T MRI system (GE Discovery 750w, General Electric, Boston, MA, USA). The protocol included a high-resolution structural scan (Sag 3D T1) with slice thickness of 1 mm, matrix 256 × 256, TR (relaxation time) 7.2 s, TE (echo time) 2.3 s, and flip angle 12 • , FOV 24, 368 slices and resting-state functional scan-2D echo-planar imaging (EPI), with slice thickness 3 mm, matrix 64 × 64, TR 2000 ms, TE 30 ms, 36 slices, flip angle 90 • , FOV 24, a total of 192 volumes. Before the EPI sequence, subjects were instructed to remain as still as possible with their eyes closed and not to think of anything in particular.

Connectivity Matrix Calculation
The normalized functional MRI volumes extracted with the help of SPM12 and the "spm_data_read" function were parcellated into 166 regions according to the automated anatomical labeling atlas AAL3 [29]. We have chosen the AAL atlas because it is the most commonly used parcellation scheme in functional network studies [30]. To estimate the connectivity between the regions of interest, we calculated an average BOLD time series x i (t) (across the voxels in each parcellation i) and Pearson correlation coefficients for all pairs of the mean parcellation activities. The Pearson correlation coefficient measures the linear relationship between two random variables and is good for low-frequency processes, which also include fMRI [31]. Connectivity Matrix calculation from the averaged activity time series was performed with the help of Matlab statistics "corrcoef" function. Thus, we obtained for each subject a 166 × 166 symmetric connectivity matrix R. Each cell of the connectivity matrix (r i,j ) represents the strength of the connection (or edge) between two parcels: Here, n is the length of the x time series, and x is the mean of the x time series.

Network Measure
To characterize the identified brain networks with meaningful and simplistic computable indicators, we used the following measures implemented in the Brain Connectivity Toolbox [32,33]: (1) The weighted undirected clustering coefficient [34], which is typically used as a measure of the prevalence of node clusters in a network. For the particular node i local undirected clustering coefficient of the network can be calculated as follows: where l i = 166 ∑ j=1 r i,j is the total weight of the relationship of the i-th node, 166 is the number of regions. The mean value for the whole network will be: (2) The weighted undirected eigenvector centrality [33,35], which measures a node's importance while considering the importance of its neighbors. This measure is defined as the eigenvector, v, associated with the largest eigenvalue λ of the correlation matrix, and can be written as: In analogy with (3), the mean centrality for the network will be denoted as <θ>.
(3) The node strength [32], which is the measure representing the sum of the edge weights: In analogy with (3), the mean node strength for the network will be denoted as <Ns>. These measures address network generative principles, order dependencies, and community structure of networks. These measures are implemented as a diagnostic tool for explicitly linking micro-scale features of network organization to macro-scale characteristics of neurophysiological dynamics [33].

Statistical Analysis
Statistical analysis of the demographic and clinical characteristics of the participants were performed using IBM SPSS 28.0 for Windows. The significance level was set to p < 0.05. We employed the Chi-square test for the categorical variables (sex and education). We tested the normality of the distributions of the continuous variables with the Kolmogorov-Smirnov one-sample criterion, which did not verify the normality for age and MADRS scores. In that regard, we employed the two-sample Kolmogorov-Smirnov nonparametric test to assess the differences between the control and patient groups.
The false discovery rate (FDR) method [36] with t-test implemented in NBS software [37] was used to assess the differences between the obtained correlation matrices for the two groups (control vs. patients) with a significance level of p = 0.0001 and 100,000 permutations. We included diagnoses as dummy-coded covariates in the t-test to control the presence of two depressive subgroups (MDD and bipolar).
For independent samples t-test was performed on the network mean values <C clust >, <θ>, <Ns>. Obtained network measures were tested for being normally distributed with the Kolmogorov-Smirnov one-sample criterion.
We applied permutation-based statistical testing [38] to test the significance difference (control vs. patients) between the distributions of the network measures by nodes. We revealed the node clusters, which were significantly different between healthy controls and depression groups based on their local values of C clust,i , θ i , Ns i . Order in the procedure is defined by the neighborhood of the anatomical parcellations, this allows us to analyze the spatial structures of the obtained clusters. For the visualization of the network structures, we used the BrainNet Viewer [39].

Linear Discriminant Analysis
To evaluate the diagnostic value of the identified differences in functional networks between the groups, we applied a linear discriminant analysis (LDA) (using the Classification Learner toolbox, Matlab) to data feature vectors from the control and depression groups. Linear discriminant analysis is widely used for diagnostic purposes; the multivariate supervised classification method sorts objects of the study into groups by finding linear combinations of a number of variables [40]. The tested variants of feature vectors included: full functional connectivity (FC) matrices, FC matrices with only connections that were considered significant through the FDR analysis, node clustering coefficients, node strength, and node centrality. Statistical significance of the results was ensured by the use of the nested k-fold (k = 10) cross-validation with each test run 1000 times. Names of classes used in the algorithm for training and classification were specified as a categorical set including "control" and "depression". Cost of misclassification of a point was set to the default Cost(i,j) = 1 if i = j, and Cost(i,j) = 0 if i = j. A linear coefficient threshold was set to zero. The discriminant type was set, recommended by the toolbox "pseudolinear" type to avoid problems with zeros and negative values in the predictors set. Other parameters including the enforced amount of regularization or prior probabilities were not applied. The chance level of the LDA for the considered problem is 50%.

Connectivity Analysis
Our analysis identified a total of 2673 connections and 152 nodes, which were significantly stronger in the control group compared to the depression group. We present the top 30 of those connections in Table 2 and in Figure 1.
Permutation-based statistical testing revealed 5 positive and 2 negative clusters of nodes in the distribution of node eigenvector centrality; 4 positive and 1 negative clusters in the distribution of node strength and 2 positive and 1 negative clusters in the distribution of clustering coefficient (see . A positive cluster corresponds to a set of neighboring nodes in which the corresponding network measure is significantly higher in the control group, while for a negative cluster, the network measure is significantly higher in the MDE group. To construct specific networks, we considered only nodes included in the significant clusters (separately for positive and negative ones) and only significant connections between such nodes that were found in the previous step with the FDR method. Specific networks for the significant positive clusters for all considered network measures are presented in a unified manner in Figure 2, and the specific networks for negative clusters are shown in Figure 3; Tables A1 and A2 in the Appendix A contain information about all connections belonging to these networks. The first two positive clusters are very similar for all network measures considered, the 2nd and 3rd are for eigenvector centrality and node strength, and the 5th is unique for eigenvector centrality. The first negative cluster is similar for all network measures considered and the 2nd is unique for eigenvector centrality.  Moreover, Clusters 2-4 are common for the node centrality and node strength in cators.  Table 3). Yellow clusters are similar for all considered network measures; green clusters similar for eigenvector centrality and node strength; blue clusters are unique for eigenvector ce trality.   Table 3). Yellow clusters are similar for all considered network measures; green clusters are similar for eigenvector centrality and node strength; blue clusters are unique for eigenvector centrality.
Moreover, Clusters 2-4 are common for the node centrality and node strength indicators.  Table 3). Yellow clusters are similar for all considered network measures; green clusters are similar for eigenvector centrality and node strength; blue clusters are unique for eigenvector centrality.   Table 3). Yellow clusters are similar for all considered network measures and blue clusters are unique for eigenvector centrality.
For the comparison Control > MDE, the main hubs common for the considered network measures are the lingual gyrus, the superior occipital gyrus, and the middle occipital gyrus.
Moreover, Clusters 2-4 are common for the node centrality and node strength indicators.

Linear Discriminant Analysis (LDA)
Employing LDA, we explored the accuracy of the diagnostic classification based on different connectivity measures. As seen in Table 6, all methods were suitable for classification purposes (precision of at least 80%), although the best results obtained from the LDA were full functional connectivity matrices. An illustration of the results are given in Figure 4, where the ROC curves are presented. ters (see Table 3). Yellow clusters are similar for all considered network measures and blue clusters are unique for eigenvector centrality.

Linear Discriminant Analysis (LDA)
Employing LDA, we explored the accuracy of the diagnostic classification based on different connectivity measures. As seen in Table 6, all methods were suitable for classification purposes (precision of at least 80%), although the best results obtained from the LDA were full functional connectivity matrices. An illustration of the results are given in Figure 4, where the ROC curves are presented.

Discussion
Our current study resulted in the following major findings: (1) patients with depression demonstrated decreased functional connectivity within as well as between different brain regions such as precuneus (PreCu), cuneus (Cu), superior occipital gyrus (SOG), lingual gyrus (LG), fusiform gyrus (FG), cerebellum, along with limbic structures including the hippocampus (Hipp) and cingulate gyrus; (2) two positive clusters (with higher measures in HC as compared to MDE patients) were common for all three network measures including node eigenvector centrality, node strength and clustering coefficient

Discussion
Our current study resulted in the following major findings: (1) patients with depression demonstrated decreased functional connectivity within as well as between different brain regions such as precuneus (PreCu), cuneus (Cu), superior occipital gyrus (SOG), lingual gyrus (LG), fusiform gyrus (FG), cerebellum, along with limbic structures including the hippocampus (Hipp) and cingulate gyrus; (2) two positive clusters (with higher measures in HC as compared to MDE patients) were common for all three network measures including node eigenvector centrality, node strength and clustering coefficient (the first cluster included mainly occipital brain regions-Cu, LG, middle and SOG while the second encompassed parts of the vermis); (3) another two positive clusters were common for both centrality and strength measures (the first one-middle cingulum bilaterally and left posterior cingulum, the second one-right anteroventral and bilateral lateroposterior thalamus); (4) one negative cluster (higher in MDE group) encompassing mainly orbitofrontal regions was common for all network characteristics while another one (with thalamic nodes) was featured solely by eigenvector centrality; (5) the LDA demonstrated that the full-connectivity matrices had the highest precision in differentiating between depression and health whilst the clustering coefficient was the least suitable measure. The significance of these findings will be discussed in the following lines.
The most significantly different connection was the one between the left and right PreCu. The function of the precuneus at rest is traditionally linked to the default mode network (DMN), which is responsible for internally oriented attention and self-reference [41]. In this context, the current findings of changed connectivity within the PreCu is not surprising and is seen in many resting-state studies of depression where the DMN, as well as its subregions, is found to be hyperactivated and hyperconnected [42,43]. This is usually linked to the clinical features of depression as a state of increased internalization (including ruminations) [44]. Notably, the direction of the change in our sample is opposite (hypoconnectivity between left and right PreCu) which is in line with the most recent meta-analytic studies finding reduced DMN connectivity especially in recurrent depression as is the case of our patient group [45].
We found no significant difference in the global (network-wide averaged) eigenvector centrality, node strength, and clustering coefficient measures between the control and MDE groups. Thus, in terms of global network topology, the characteristic functional networks do not differ between the groups under consideration. The main differences are observed at the local level, as indicated by the significant clusters in the network measure distributions (see Tables 3-5).
The first two positive clusters (Control > MDE) in the network measure distributions, in which eigenvector centrality, node strength, and clustering coefficient are higher in the control group, include bilateral lingual, superior, and middle occipital gyri and cerebellar regions (vermis). In terms of the network measures, a higher local clustering coefficient means that short-range (local) connections prevail over long-range (global) ones in these brain areas in the control group, and network clusters are formed there. A higher eigenvector centrality and node strength in the same nodes in the control group corresponds to a stronger integration of emerging clusters in the network and stronger connections of these nodes (hubs) with other large hubs. The remaining three positive clusters (the cingulate, thalamic nodes and the inferior frontal gyrus) are characterized by increased eigenvector centrality and node strength (only for the 3rd and 4th clusters) in the control group. This indicates that these nodes form larger and more strong integrative network hubs compared with the MDE group. These conclusions are also supported qualitatively in Figure 2.
Of note, all considered measures pointed to the role of the left LG as a major hub in the occipital cluster, with a number of connections demonstrating a significant difference in the depressed as compared to healthy individuals. In accordance with these results, depression has been linked to impaired static and dynamic FC of the LG. The role of the left LG in depression has been reported in a most recent functional connectivity study assessing the effect of childhood trauma [46]. The authors report that the FC changes of the left LG were not affected by the presence or absence of traumatic events, which may reflect a general vulnerability to depression. In support of this notion is a recent study by Wang et al. on electroconvulsive treatment of depression. The authors found changes in functional connectivity of the lingual gyrus to be persistent before as well as after the procedures [47]. The other common hub for all three network measures in our study was the vermis. Apart from the well-known motor functions, there is growing evidence for the involvement of different parts of the vermis in higher order functions including cognitive and emotional processing [48]. Although long neglected, the role of the cerebellum in emotion has been now reestablished also in a recent consensus paper [49].
The links with psychiatric disorders are not surprising. Cerebellar gray matter reductions as well as decreased activity and connectivity of this region have been reported in depression [25,50]. Some authors proposed that impaired cerebellar function contributes to abnormalities in predictive coding and homeostatic dysregulation in depressive disorder [51].
The third cluster which demonstrated significantly different distribution of node strength and node centrality encompassed the mid-cingulate cortex(MCC) which is divided into the posterior MCC involved in multisensory orientation of the head and body in space while the anterior MCC is involved in action-reinforcement associations and selection based on the amount of reward or aversive properties of a potential movement. MCC contributes to cognitive control and decision making [52]. The anterior subregion also has high dopaminergic afferents and high dopamine-1 receptor binding and is engaged in reward processes [53]. Emotional n-back tasks elicited differential activation of the posterior MCC in unipolar compared to bipolar depression [54]. In addition, MDD patients as compared to healthy individuals failed to activate MCC when an emotional stimulus was paired by a neutral one [55].
The fourth positive cluster yielded by the comparison HC vs. MDE patients included different areas of the thalamus with bilateral lateroposterior and mostly right-sided anteroventral engagement. The anteroventral/anteromedial together with the mediodorsal nuclei play important roles in connecting subcortical limbic structures (amygdala) to the limbic cortex (anterior cingulate and orbitofrontal cortex) [56]. In addition, the cerebellothalamo-cortical loops are implicated in emotion regulation and subjective sense of control [57]; and aberrant intrinsic FC of the thalamocortical pathway was associated with depression [58].
The first negative cluster (MDE > Control), which is the same for all considered measures, includes the following brain areas: superior frontal gyrus (medial orbital), rectus gyrus, and medial and anterior orbital gyri. Thus, in the MDE group, more developed network clusters are formed near these nodes, while they are also more strongly integrated into the whole network, being large hubs (see also Figure 3). The second negative cluster is unique for the eigenvector centrality measure and includes thalamic nodes. Note that a significantly large number of positive clusters in comparing the control vs. MDE groups is a characteristic feature of the functional network.
The orbitofrontal cortex (OFC) plays a key role in emotion, by representing the reward value of the goals for action [59]. Its involvement in depression has been supported by a number of structural and functional imaging studies [10,60]. The ENIGMA consortium found that MDD patients had thinner gray matter in OFC [61]. Moreover, one of the most recent hypotheses about the development of depression, namely the non-reward attractor theory of depression, assigns a major role of the lateral OFC in the pathophysiology of the disorder [62]. According to the author the lateral orbitofrontal cortex non-reward system triggers negative cognitive states, which in turn have positive feedback top-down effects on the orbitofrontal cortex non-reward system. Our results of increased network measures of the anterior OFC (corresponding to the functional division of lateral OFC] hub in depressed patients is in line with this recent theory.
On the other hand, we know that the medial OFC is activated by rewarding and subjectively pleasant stimuli and MDD has been found to be characterized by reduced functional connectivity of the mOFC with the parahippocampal gyrus, which is in line with the clinical manifestations of anhedonia (reduced anticipatory and consummatory pleasure [63]. Our results support previous findings of depression being characterized by disturbances of both lateral and medial part of the OFC.
The most important to the clinical reality findings of the current study was derived from the linear discriminant analysis. It demonstrated that the clustering coefficient was the most ineffective measure, while full-connectivity matrices, as well as those with only the significant connections identified in advance, were the most precise in differentiating between depression in patients and healthy individuals. These measures reached precision levels of 97% and 94%, respectively. Thus, the connectivity matrices outperformed the network-specific features of node strength and node centrality. We can speculate that in order to demonstrate a meaningful diagnostic value, the resting-state connectivity feature of depression should be considered as a whole and not reduced to separate network measures.
Earlier connectivity-based classification studies focused on specific regions or networks, some of them reaching very good accuracy levels of around 90% [64]. Later, wholebrain connectivity analysis was favored, and the prediction levels reached 94% [65]. Some of the most recent studies adopted a fusion strategy using different connectivity features (intrinsic, dynamic FC, effective connectivity) that could distinguish between MDD and controls with an accuracy of 90.91% and an AUC of 0.895 [66]. In this regard, our study yielded one of the highest performances of the classifier based on whole-brain functional connectivity analysis.

Limitations
This study has several limitations that must be admitted. Firstly, the sample size is relatively small, although comparable to other single-site studies. Secondly, the patient group is heterogeneous in terms of diagnosis (major depressive and bipolar disorder) and this may have influenced the results since both common and distinct activity and connectivity patterns have been demonstrated [67,68]. However, we have studied dysfunctional connectivity as a state dependent measure in major depressive episode. Unlike the trait (or state-independent indicators), clinical features as well as the underlying functional networks impairment on the level of the episode or the syndrome are more or less shared between BAD and MDD. Thirdly, the fact that all patients have been on stable antidepressant medication prior to inclusion might have contributed as well to our findings. Future research should address those limitations by increasing the number of subjects, exploring the two types of depression separately, and including samples of non-medicated patients.

Conclusions
In the recent decade, there is a major paradigm shift toward the network-connectivitydriven classification of mental disorders [69]. Our work contributes to this field by providing insights into impairment of resting-state network hubs elucidated by functional connectivity measures such as node centrality, node strength and clustering coefficient. This approach delineated brain structures, including the lingual gyrus, cerebellum, midcingulate cortex and thalamus, as being more prominent in healthy subjects as compared to depression where the orbitofrontal cortex emerged as a key node. In addition, the connectivity matrices proved to be suitable for differentiating between patients with MDE and HC with fairly high precision levels. Further integration of these results with clinical measures, structural and diagnostic task-related functional MRI on the level of multivariate analysis may outline new approaches and criteria to the definition of mental disorders [70][71][72].  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
All data is available from the authors upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest. Table A2. List of significant (found with the FDR method) between negative clusters nodes (See Figure 3) for the node centrality measure.