A Resting-State Brain Functional Network Study in MDD Based on Minimum Spanning Tree Analysis and the Hierarchical Clustering

1School of Information Science & Engineering, Lanzhou University, Lanzhou, China 2International WIC Institute, Beijing University of Technology, Beijing, China 3TheThird People’s Hospital of Tianshui City, Tianshui, China 4Lanzhou University Second Hospital, Lanzhou, China 5Beijing Anding Hospital of Capital Medical University, Beijing, China 6Computer Systems Institute, ETH Zürich, Zürich, Switzerland


Introduction
Major depressive disorder is a global mental disorder and has an unfavourable influence on physical and psychological health [1].In addition to profound personal suffering, MDD patients lack the necessary social and occupational functioning [2].Moreover, The World Health Organization predicted that depression would become the second leading cause of illness by the year 2020 [3].In this light, exploring the neurobiological signature of MDD from multiple imaging modalities was considered to sharpen the reach of depression and develop treatments, including electroencephalogram (EEG), magnetoencephalogram (MEG), functional magnetic resonance imaging (fMRI), positron emission tomography (PET), and single photon emission computed tomography (SPECT) [4].In recent years, the research results of MDD based on different approaches had been presented substantially such as frontal EEG asymmetry, "small-word" network characteristics, and increased/disrupted cognition connectivity network [5][6][7][8].These results revealed neurophysiology characteristics in different aspects for depression disease and made a great contribution to the study of the depression.However, there were disputes and contradictions in these results due to the differences of subjects, experimental environment, methods, and other restrictions.So, more methods and techniques are expected for exploring MDD.

Complexity
The human brain is a complex system, characterized by its dynamical neural communications and mutual interactions based on synchronous oscillations among different brain areas [9,10].In the human brain, oscillatory patterns reflect the activity of the brain and provide reliable markers of the brain function or dysfunction [11,12].Therefore it is a direct and effective way to explore the brain activity based on the brain oscillatory nature.In present imaging modalities, electroencephalogram (EEG) is regarded as a convenient technology to reflect the comprehensive electrophysiological activity of neuron populations [13].In recent years, it has been widely used in biomedical fields to estimate the oscillatory patterns on account of higher temporal resolution and noninvasiveness in mental disorders such as major depressive disorder (MDD), Alzheimer's disease, and schizophrenia disease [14,15].EEG oscillations are rhythmic electrical events coming from the brain and can be used to define the interaction of different brain regions [16].Because of this feature, it is suggested that the information processing of the brain can be reflected in characteristic EEG oscillation rhythms [17].Using this approach, a large number of findings were presented in the study of the depression.Some results based on EEG oscillations demonstrated that patients with MDD had more frontal theta, alpha, and beta oscillations [18][19][20].Using EEG oscillations, Lee et al. [21] suggested that the brain affected by a major depressive disorder showed a slower decay of the long-range temporal (auto)correlations (LRTC).Recently, an EEG oscillations study on MDD reported that depressive brain was manifested in the superposition of distributed multiple oscillations [22].
At present, studying functional interactions between brain regions plays a vital role in understanding the dynamic interactions between the neural systems [23].For defining the interaction of different brain regions, the synchronization of EEG oscillations is an important and effective indicator which can be estimated by EEG coherence [24].And EEG coherence is conceptualised as the correlation in the time domain between two signals in a given frequency band [25,26].The high EEG coherence reflects synchronized neuronal oscillations between different brain areas, whereas low EEG coherence represents independent active neuronal [25].This approach has been applied for evaluation and auxiliary diagnosis of various mental disorders, including Attention-Deficit/Hyperactivity Disorder, Autism Spectrum Disorder, MDD, and Alzheimer's disease.In one study of MDD, EEG coherence was used to estimate the sleep EEG rhythms, which suggested that low temporal coherence in depression reflects a breakdown in the organization of sleep EEG rhythms within and between two hemispheres [27].Li et al. [7] found that the global EEG coherence of patients with MDD was significantly higher than that of healthy controls in both gamma bands.Prior EEG coherence based on discriminant function analysis (DFA) rules was used to explore possible neurophysiological differences between Asperger's Syndrome (ASP) and the Autism Spectrum Disorders (ASD) and successfully distinguished ASP and ASD populations [28].Using EEG source-based coherence in Alzheimer's disease (AD) showed increased delta coherences between the bilateral precentral, left supplementary motor area (SMA), and right precentral [29].
In recent years, with the analysis of the connectivity based on EEG data, graph theoretical analysis has been widely applied.With the rising interest in graph theoretical studies of brain networks, the problem of the determination of the connectivity structure in the brain has been a subject of the intense research.The approaches standing behind the connectivity have a crucial impact on the results of the study and complicate a comparison of results across different studies and different brain imaging techniques [30].Conventionally, a threshold, or a range of thresholds, is used to confirm whether connections exist or do not exist.Currently, this method is most widely used in the study of the brain functional networks [31][32][33][34].But it might be biased by the choice of the threshold due to the number of links being decided by the size of the threshold values.Importantly, most network characteristics depending on the number of links in the network would be also influenced such as clustering coefficient, characteristic path length, and node degree.It may be an assignable cause that there were some controversial findings in the study of the brain functional network.In order to avoid this bias, some of the studies used the weighted brain functional network to avoid the choice of threshold [35][36][37].Although the weighted network overcomes the problem of this subjective factor and provides a more realistic representation of functional networks, there are also problems in the weighted network.Spurious weak connections are also taken into account, potentially influencing the brain functional network.
In this context, using minimum spanning tree (MST) to represent brain networks may be one promising unbiased solution to this problem [38].The MST is a subgraph without forming cycles that connects all the nodes in the original weighted network [39].In this way, MST will obtain the same number of nodes and links, therefore enabling the direct comparison of network properties between groups and avoiding the aforementioned methodological biases and defect.Due to this advantage, the MST has been wildly used in a variety of mental diseases.Stam et al. [30] illustrated MST characterization allowed the representation of the observed brain networks and may simplify the construction of simple generative models of normal and abnormal brain network organizations.MST appeared in a variety of studies.MST was used as an elegant and sensitive method to capture subtle developmental organization changes in the brain networks of children [40].In a study of Multiple Sclerosis, findings indicate that MST network analyses were able to detect network changes in the Multiple Sclerosis (MS) patients [41].In the study of Alzheimer's disease, MST was regarded as an effective method for analyzing cortical networks [42].
The clustering and community structures have been regarded as one of the most significant features of complex networks [43].In the brain networks, clustering or community structure was defined as a subset of highly interconnected nodes which had similar characteristics [44].Moreover, we have known that brain networks demonstrate the property of the hierarchical modularity.Constructing the brain hierarchical modularity, the hierarchical clustering is a special method which characterizes major building blocks or the hierarchical modularity of brain networks, corresponding to specialized brain functions [45].Among the multitudinous hierarchical clustering algorithms, a tree agglomerative hierarchical clustering (TAHC) method can successfully detect clusters in both artificial trees and the MSTs of weighted social networks.This method was raised in [43].Moreover, the hierarchical clustering has regarded the MST as a foundation of the complex networks [30].And the hierarchical clustering in combination with the MST has been used in different cognitive domains.Hierarchical clustering analysis of the MST showed that the connections between the parahippocampal gyrus and posterior cingulate gyrus were disrupted in Alzheimer patients [42].Impaired communication between functional clusters in AD was found in one MEG study [46].Alexander-Bloch et al. [47] found disrupted modularity and local connectivity of the brain functional networks in childhood-onset schizophrenia using hierarchical clustering analysis of the MST.
This study used eyes-closed resting-state EEG recordings involving the patients with MDD and the healthy subjects.The aim of this study was using an unbiased method to construct the brain network and then explore the network characters and clustering of nodes for both MDD and healthy groups.Although the study of "resting-state" is rather exploratory, we believe that the approach used in this study is feasible.The choice of threshold has a great impact on the properties of the brain functional network.So we used an unbiased MST method to build the main network of the brain.As far as we know, this method was first used for the depression disease.Then we estimated the leaf fraction, mean link weight, and node degree fraction of MST.The hierarchical partitioning of the brain functional network in resting-state for MDD and healthy groups was not clear.So, in this study, TAHC method was used to characterize major building blocks and hierarchical partitioning of the brain networks.In the end, we discriminated the differences between the MDD group and the healthy group and then summarized and discussed our findings.

Subjects.
The study was approved by the local ethics committee.Written informed consent was obtained before the study began.Twenty-three patients with MDD (13 males and 10 females, right-handed) were recruited from the Lanzhou University Second Hospital.The mean age of the MDD group was 33.17 ± 19.83 years.All patients had no history of the manic episode.14 healthy subjects (7 males and 7 females, right-handed) were recruited from the society with the mean age of 31.29 ± 21.71.To ensure the effectiveness of this study, the participants were aged between 18 and 55.Besides, primary or higher education level was required.Strict exclusion criteria were enforced before the experiment, and exclusion criteria for all participants included past history or the presence of any medical or neurological disorders, presence of drug or alcohol abuse, and past head trauma with loss of consciousness.Before the experiment, all participants participated in an interview in which the Mini and PHQ-9 [48] were administered with the help of an experienced clinical psychiatrist.Mini was used to ensure the correctness of the classification.The score of PHQ-9 was used to evaluate state anxiety, general anxiety, and depression levels.The mean PHQ-9 score of MDD group was 17.10, and the healthy group was 2.57 ( = 11.504, < 0.001).
In addition, all participants gave informed consent and were rewarded for their participation.

EEG Data Processing.
The experiment was conducted in a quiet and dim light room kept away from electromagnetic interference.Participants were required to have a seat on a wooden chair comfortably with eyes closed but keeping awake during the EEG recording.They were also asked to avoid blinking and making movements.

Coherence.
The coherence is defined as the spectral cross-correlation between two signals normalized by their power spectra [7].There are different measuring methods that analyze the coherence from different pairs of electrodes per frequency.In this study, the magnitude-squared coherence (MSC) was calculated for a particular frequency  between two given EEG signals  and .
() is the power spectral density (PSD) estimate of  at the frequency of  and   () is the cross PSD estimate of  and  at the frequency of , using Welch's averaged, modified periodogram method.The value of the MSC ranges between 0 and 1, where 0 represents no coherence and 1 indicates maximum linear interdependence between two signals.
In this study, we calculated the coherence between each possible pair of 72 EEG channels with the respect to each single frequency.A square 72 * 72 coherence matrix was obtained for each participants (72 was the number of the chosen EEG channels), and each element in the coherence matrix indicated the coherence between the corresponding two electrodes.The MDD coherence matrix was defined as the mean of all coherence matrices of the patients with MDD, and the definition of coherence matrix for healthy group was the mean of all coherence matrices of the healthy participants.The global coherence was defined as the mean value of all elements in the coherence matrix.It indicated the level of interdependence of the whole brain.In order to confirm which frequency bands have significant difference rapidly, we calculated the global coherence between two groups in each single frequency within 3-35 Hz.

MST.
The minimum spanning tree is a simple acyclic connected subgraph of the original weighted network that can be used to direct comparison of networks with the same number of nodes and simplifies the network characterization.Although it is not a priori given, it will still capture most of the important topological information in the original network [30].We constructed the MST based on the aforementioned EEG coherence matrix by employing Kruskal's algorithm [49].The details of this algorithm used in this study were as follows: (1) ordering the elements of the EEG coherence matrix in a degressive order; (2) linking the  nodes with maximal EEG coherence until all the nodes being linked in a loopless subgraph consisting  − 1 edges; (3) skipping the link, if adding this link leading a circle.
In this study, we used 72 channels.So the number of the nodes in the topology of MST was 72 and the number of edges was 71.We also analyzed the topology properties of the MSTs for the two groups.Leaf fraction and the mean weight of all links included in the MST were calculated to evaluation of the MST topology.After that, statistical analyses were used to assess the credibility of differences.

Hierarchical Clustering.
The clustering is an effective method to explore nontrivial information in the network.MSTs contain most of the information about the underlying clusters of the original weighted networks.In this study, TAHC method was used for the detection of the clusters in the aforementioned MSTs.The summary of TAHC algorithm [43] is as follows: first we use geodesic distances between all possible pairs of nodes of the given graph as an input to the agglomerative hierarchical clustering algorithm.Next, compute the similarity between every node pairs in an MST based on geodesic distances.Then find the most similar pair of clusters and merge them into a single cluster.Finally, recalculate similarities between the new cluster and each of the old clusters based on average-linkage clustering and remerge clusters until all nodes are merged into a single cluster.
The geodesic distance between two nodes in a tree is equal to the number of links in the shortest path.So, geodesic distances were calculated using Dijkstra Shortest Path algorithm [50].Geodesic distances between all possible pairs of nodes in a graph constituted the geodesic distance matrix  which was a weighted matrix.In the geodesic distance matrix , each node corresponds to a row vector.Based on , we calculated

Global coherence
Figure 1: The coherence-frequency graph of MDD group and healthy group.The horizontal axis denoted the frequency and the vertical axis denoted the global coherence.There was the significant difference in the theta band (4-8 Hz).In alpha (8-13 Hz) and beta bands (13-30 Hz), there was no diacritic difference in both groups.
vector similarities using Spearman's rank correlation between all row pairs of .
is Spearman's rank correlation, and the value of the MSC ranges between 0 and 1.The larger the value is, the higher the similarity of the two nodes is. is the number of the entries in a row.  is the difference between the two rows of each observation.
In our study, we analyzed the hierarchical clustering organization of the group-level MST for each group in interested frequency bands.Thus, each group-level dendrogram obtained by the TAHC method corresponds to each grouplevel MST.The distribution of nodes in the hierarchical clustering was described based on a global electrode graph.

Global Coherence.
A large number of researches indicated that most of the EEG signals were at low frequencies when the brain was in the resting-state.In order to explore which frequency bands contain significant differences between the MDD group and the healthy group rapidly, we calculated global coherence in each single frequency within range of 3-35 Hz and used a frequency-coherence line chart to describe the relationship between the global coherence and the frequency, as seen in Figure 1.In this graph, we can see that the higher coherence appeared in the alpha band (8-13 Hz) for both groups.And the global coherence of the MDD group is significantly higher than that of the healthy group in the theta band (4-8 Hz).In alpha (8-13 Hz) and beta bands (13-30 Hz), there is no diacritic difference in both groups.

Difference in Coherence.
Because the difference obtained from the global coherence is not obvious in the alpha and beta bands, we paid more attention to the theta band.We calculated the coherence of each pair of the channels.The coherence matrices were plotted in Figure 2. The coherence matrices of both groups showed a complex but rather similar pattern, with various regions of high and low levels of the interdependence.There was no evident difference of the distribution for highlight areas in both groups.Compared with the healthy controls, the patients with MDD had more strengthened coherence values in some regions corresponding to the red areas.In order to explore the distribution of this difference, we obtained a difference matrix using the coherence matrix of the MDD group to subtract the coherence matrix of the healthy group.We plotted this matrix in a 3D graph as seen in Figure 3.The topographic maps were plotted by the threshold of 0.12 and 0.14, which could clearly indicate the difference between two groups in the theta band according to our above results.At the threshold of 0.12, most of the links were distributed in the left hemisphere of the brain except the front regions.There were also some links in the right-temporal region.At the threshold of 0.14, most of the links were distributed in the left hemisphere of the brain especially in the parietal and temporal regions.There were few links in the right hemisphere of the brain.
Coherence analysis results showed that the depressed group had significantly higher coherence in the left hemisphere of the brain especially in parietal and temporal regions compared with the healthy controls in the theta band.

Characteristics of MST in Both Groups.
In order to avoid the biase coming from the choice of the threshold, MST was used to construct the brain functional network.In one extreme, all nodes are connected to two other nodes, with the exception of the two nodes at either end, which have  only one link.The other extreme is a star.In this star, there is one central node to which all other nodes are connected with one link.This tree configuration would have more leaf nodes.Between the two extremes of one path and one star, many different types of the tree configuration are possible.We calculated the leaf fraction for both groups.The leaf fraction of the MDD group was 0.5984 and for the healthy group was 0.5437 ( = 3.1368,  = 0.0035).Compared with the healthy group, MST of the MDD group tended to a star-network which had more leaf nodes.We also calculated the mean weight of all links.The mean weight of the MDD group was 0.8464 and for the healthy group was 0.7967 ( = 2.2168,  = 0.0335).Then the degree fraction of nodes was presented in the brain topographic mappings.Figure 4 showed the results.
In the topographic map, we could see that clearly the degree fraction of the MDD group in the temporal regions in two hemispheres of the brain was higher than that of the healthy group.But the number of center nodes in the MDD group was less than that of the healthy group.

Hierarchical Clustering Analysis
. TAHC algorithm was used to construct the hierarchical clustering of the brain functional network for both groups.Figure 5 showed the results.In the hierarchical clustering of the MDD group, it tended to cluster according to the physical structure of the brain, and the clusters corresponded to the left and right hemispheres.In the healthy group, the clustering tended to the functional structure of the brain, corresponding to the different functional areas.In order to explore the clusters detailedly, we plotted the distribution graph of the clusters in Figure 6.The hierarchical clustering was integrated into 6 clusters for the MDD group and the healthy group.We found a significant difference in the front region.In this region, there was a cluster which contained a large number of nodes in the distribution graph of the healthy group.However the MDD group had not.

Discussion
The results showed an increased global coherence of the MDD controls compared with the healthy controls in the theta band.In the left hemisphere of the brain, the MDD group had higher coherence, especially in the parietal and temporal regions.Several studies have reported that there was abnormal EEG absolute power, hemispheric asymmetry, or coherence with depression patients in some specific frequencies [51].Hinrikus [52] found that the controls with MDD had the increased coherence between some brain regions.Tucker and Dawson [53] also reported that the parietal and temporal cortical regions were more activated in the left hemisphere of the depressive patients based on the coherence measure.
In other frequency ranges, the increased coherence was also found.Fingelkurts et al. [54] reported the increased synchronization was observed in the EEG alpha and theta bands in the patients with MDD.And an increased topographic EEG coherence in the frontal brain areas in the MDD group in the EEG alpha, beta, and theta bands was reported by Leuchter's team [55].However, the findings in the previous studies are inconsistent and sometimes even contradictory to each other.Knott et al. [19] found the decreased coherence in MDD subjects compared to normal controls.It may be due to the limited number of individual electrodes or the improper measurement method.For the increased global coherence of the MDD group, some researchers interpreted it as adaptive and compensatory mechanisms aimed to overcome the deficient semantic integration [55].Our findings further support this adaptive and compensatory mechanism.However, it is not clear why the higher coherence was obtained in the left hemisphere of the brain.But the hemispheric asymmetry in the emotion processing has been observed, and the functional complementation of left and right hemispheres is important for adaptive emotion regulation [56].Deficit of emotion regulation ability is the main manifestation of the patients with MDD [57].So we guess it may be related to the dysfunction of the left hemisphere.Future studies about function and activation of the left and right hemispheres are needed to reach more definitive conclusions.The functional network of the patients with MDD has been widely explored and studied.And a large number of great achievements have already been obtained.But, due to the differences in environments, subjects, and methods during the experiment and data analysis process, different  or even opposite results were obtained in the research of the MDD brain functional networks.Some teams discovered increased brain functional connectivity in the patients with MDD [54,[58][59][60], and decreased brain functional connection was also found in Yang et al. [61] and Wang et al. [62] teams.For the functional connection, the choice of the threshold was commonly used to convert the coherence into functional connection.It is a reason that could not be ignored for emergence of the different conclusions.In this study, we used an unbiased method MST to construct the main brain functional connection.This method has been used in multiple mental disorders.As far as we know this was the first time for this method to be used for the depression disease.In a previous study, it suggested that more random networks showed low clustering and a short path length, corresponding to MSTs' shorter diameters and higher leaf numbers [40].Our finding that leaf fraction of the MDD group was higher than the healthy group indicated a shift toward randomization in the brain networks of the MDD group.Similar conclusions have been mentioned in the study of brain functional networks.In sleep neuronal functional networks of depressed patients, Leistedt et al. [63] indicated the functional reorganization of depressed patients lost "small-world network" (SWN) characteristics.Zhang et al. [64] and Li et al. [7] teams also indicated the MDD patients showed a shift toward randomization in MDD brain networks compared with the healthy controls.For the MST, the higher the degree of a node is, the more important this node is during the brain information processing.Our finding about degree fraction indicated that temporal regions played an important role in information processing of the MDD group.Previous studies had the analogous conclusion.An EEG source location study suggested that larger degree fraction in temporal regions of the MDD subjects may be related to the dysregulated temporal pole activity [14].The prefrontal cortex mediates the control of high-level cognitive functions and is associated with the regulation of many aspects of the affective system [65].This idea had been supported in neuroimaging studies.For the study of depression, the depression is known to involve a disturbance of mood and impaired cognitive functions [66,67].In previous studies, the frontal regions obtained more attention and exploration and played an important role in the development of the depression [68].So far, many results reveal that major depression can be distinguished by specific histopathology of both neurons and glial cells in the prefrontal cortex [69].Biver et al. [70] found frontal metabolic disturbances in the unipolar depression.Disruption of paralimbic pathways linking frontal cortex in secondary depression was indicated in [71].But the most mentioned finding was the frontal brain asymmetry [72][73][74][75], which was described with greater activation in the right compared to the left frontal lobes [76].The asymmetric frontal cortical activity in the MDD group had been widely presented not only in alpha band but also in the theta band [24,77].Besides EEG study, a combined MEG, PET, and rTMS Study also pointed out that prefrontal left-right functional imbalance and disrupted prefrontothalamic circuitry were plausible mechanisms for the depression [78].Hierarchical clustering was a useful method to divide the brain functional network into several submodules.The nodes in the same submodule had a strong similarity, and the information processing among these nodes was very efficient.Importantly, the submodules typically corresponded to functional systems of the brain.In our study, we found that there was a cluster which contained a large number of nodes in the healthy group.However the MDD group had not, and the nodes in frontal regions of the MDD group were divided into two clusters in the left forehead and the right forehead, respectively.It was the performance of the forehead imbalance in patients with depression.This result indicated that there was a left-right functional imbalance in the frontal regions for MDD controls.

Conclusion
In conclusion, abnormally increased EEG coherence of the MDD group was found in the theta band, and the higher coherence was described in the left hemisphere of the brain especially in the parietal and temporal regions.An unbiased method of MST was used to construct the brain functional networks for the MDD group and the healthy group.The higher leaf fraction and mean weight were found in the MDD group.This finding indicated a shift toward randomization in the brain networks of the MDD group.Additionally, the hierarchical clustering opened up a new way for obtaining the characteristics of the brain functional network.The results that the MDD controls lose a frontal clustering indicated that the MDD group lacked the coordination in the forehead.A possible disadvantage of the MST approach is that it may miss the information about the network topology and it may contain the weaker connections in this brain functional network.Then, in the further study, we will try to use two levels' MST to construct the brain functional networks to overcome the above problems.

Figure 2 :Figure 3 :
Figure 2: The coherence of MDD group and healthy group in the theta band.The size of the coherence matrix was 72 * 72.In the matrix map, each chromatic point represented the coherence of two corresponding channels.The horizontal and vertical axes denoted 72 channels.
Degree fraction of the MDD groupDegree fraction of the healthy group

Figure 4 :
Figure 4: The brain topographic mappings of the degree fraction.The depth of the color represents the degree size of the nodes.

Figure 6 :
Figure 6: The distribution of the hierarchical clustering graphs under the condition of six clusters.