Intrinsic Organization of Occipital Hubs Predicts Depression: A Resting-State fNIRS Study

Dysfunctional brain networks have been found in patients with major depressive disorder (MDD). In this study, to verify this in a more straightforward way, we investigated the intrinsic organization of brain networks in MDD by leveraging the resting-state functional near-infrared spectroscopy (rs-fNIRS). Thirty-four MDD patients (24 females, 38.41 ± 13.14 years old) and thirty healthy controls (22 females, 34.43 ± 5.03 years old) underwent a 10 min rest while their brain activity was recorded via fNIRS. The results showed that MDD patients and healthy controls exhibited similar resting-state functional connectivity. Moreover, the depression group showed lower small-world Lambda (1.12 ± 0.04 vs. 1.16 ± 0.10, p = 0.04) but higher global efficiency (0.51 ± 0.03 vs. 0.48 ± 0.05, p = 0.03) than the control group. Importantly, MDD patients, as opposed to healthy controls, showed a significantly lower nodal local efficiency at the left middle occipital gyrus (0.56 ± 0.36 vs. 0.81 ± 0.20, pFDR < 0.05), which predicted the level of depression in MDD (r = 0.45, p = 0.01, R2 = 0.15). In sum, we found a more integrated brain network in MDD patients with a lower nodal local efficiency at the occipital hub, which could predict depressive symptoms.


Introduction
Major depressive disorder (MDD) is a disabling disease associated with profound functional impairment [1]. According to the World Health Organization, before the COVID-19 pandemic, approximately 350 million people worldwide were suffering from depression [2]. Estimates put the rise in both anxiety and depressive disorders at more than 25% during the first year of the pandemic [3]. In addition to the core emotional symptom, cognitive subdomains, such as learning and memory, executive functioning, processing speed, and attention and concentration, are significantly impaired in patients with MDD [4]. However, due to the lack of specific diagnostic methods for depression and its various manifestations, accurate and reliable diagnosis of MDD mainly depended on the subjective assessment and clinical experience of the clinicians. Thus, the accuracy of reaching a diagnosis, relying on past experience remains debated [5]. Biomarkers offer a conceivable target for assisting in diagnosis and identifying predictors of response to various interventions [6]. Brain imaging techniques are powerful tools for tracking biomarkers, and thus have been used as an adjunct to help clinicians diagnose MDD in recent years.
Functional near-infrared spectroscopy (fNIRS) is a non-invasive, functional neuroimaging technique that could assess cerebral function by quantifying cerebral oxygenation [7][8][9][10]. This type of spectroscopy, fNIRS, penetrates organic tissues by using light sources with a spectral window of 650 to 1000 nm. The variance in absorbance is then used to compute

Hypotheses
The first hypothesis is that the dysfunction of the brain network in patients with MDD can be detected by rs-fNIRS. Thus, we sought to perform an rs-fNIRS based functional connectivity analysis to investigate the intrinsic organization of brain networks in MDD. Moreover, we have the second hypothesis that the disrupted topological architecture of the brain network can be found in patients with MDD using rs-fNIRS, which could help distinguish patients from healthy people. Therefore, we calculated graph theory-based connectivity metrics to quantitatively characterize functional connectivity in each group, expecting to find some brain indices that helps differ MDD to healthy individuals.

Participants
A total of 64 adults, including 34 participants (24 females, mean age ± SD = 38.41 ± 13.14 years old) in the depression group and 30 participants (22 females, 34.43 ± 5.03 years old) in the healthy control group, participated in the current study. All patients were inpatients who volunteered to take part in this study and provided their consent forms, and their admission dates ranged from 25 February 2021 to 31 August 2021. The patients with depression were recruited from the Hangzhou Seventh People's Hospital, while the healthy controls were recruited through advertisements. All the patients with depression met the American Psychiatric Association DSM-IV diagnostic criteria of depression, and healthy controls were interviewed using the Structured Clinical Interview for DSM-IV, nonpatient edition. Age and sex were matched between the two groups (age: t = 1.56, p = 0.12; sex: χ 2 = 2.08, p = 0.20). In the non-depressive sample, 6.67% (n = 2) of people had a doctorate degree, 23.33% (n = 7) had a master's degree, 63.33% (n = 19) had a university degree, and 6.67% (n = 2) had an associate degree. A total of 17.65% (n = 6) of the depressed sample has a university degree, 11.76% (n = 4) had an associate degree, 20.59% (n = 7) had a high school diploma, and 50% (n = 17) had a middle school degree. A total of 61.76% (n = 21) of the clinical sample had a generalized anxiety disorder; 14.71% (n = 5) of all depressive patients were diagnosed with non-organic sleep disorders and other comorbidities, such as sleep apnea syndrome (5.88%, n = 2), obsession (5.88%, n = 2) and panic attacks (5.88%, n = 2). Hamilton Depression Rating Scale (HAMD) [40] was used to assess the depression level in both groups (depression vs. control, 22.55 ± 1.27 vs. 6.66 ± 1.51, t = 17.09, p < 0.001, Cohen's d = 11.39). The study procedure was carried out following the Declaration of Helsinki and approved by the ethics committee of local hospital (No. 2021024). All participants gave their written informed consent before the experiment.

NIRS Data Acquisition
We collected our resting-state data using a multi-channel functional near-infrared spectroscopy (fNIRS) optimal imaging system (NirScan-900A, Huichuang, China), equipped with 15 detectors and 19 sources, resulting in 39 measurement channels in total to cover prefrontal, central, and posterior cortices (see Figure 1 for the locations of measurement channels). The distance of source-detector separation was approximately 3 cm. The absorption of near-infrared light at three wavelengths (730, 808, and 850 nm) was recorded with a sampling rate of 50 Hz. The light absorption data were converted into concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) using a modified Beer-Lambert law [41], with a differential pathlength factor of 6 [42]. We focused on HbO signals given the higher signal-to-noise than HbR [43]. For each participant, the rs-fNIRS data were collected for 10 min, during which the participants were instructed to keep relaxed, keep their eyes open, and remain awake. Brain Sci. 2022, 12, x FOR PEER REVIEW 4 of 16

Network Construction
Network construction encompassed two steps. First, in accordance with previous studies [47], Pearson correlation coefficients were used to quantify the relationships among time courses of every pair of channels, resulting in a 39 × 39 correlation matrix. We zeroed all the negative correlation coefficients because of their ambiguous biological meanings [48,49]. We restricted following analyses to only positive correlations. All correlation coefficients were subjected to Fisher's z-transformation to improve normality.
Second, we binarized the correlation matrix according to sparsity-based threshold [50,51]. The correlation matrix was thresholded over a range of sparsity (from 20% to 40%, with 1% step size) in order to investigate the relationship between sparsity and the network properties. Specifically, each functional connectivity matrix Zij can be converted to a binarized matrix Bij, where Bij is 1 if the value of the z value in matrix Zij is greater than a given sparsity threshold and 0 otherwise [52]. Consistent with recent recommendations [52], we also adopted a single sparsity (25%) to normalize the network metrices, thereby allowing the comparison of group differences under the same topological organizations.

Network Analysis
In the current study, we used the GRETNA toolbox [53] in MATLAB 2022a to assess the topological measures, including global network metrics (small-world Gamma, Lambda, Sigma, global efficiency, and local efficiency) and regional nodal metrics (global nodal efficiency and local nodal efficiency) for the brain network.

Network Construction
Network construction encompassed two steps. First, in accordance with previous studies [47], Pearson correlation coefficients were used to quantify the relationships among time courses of every pair of channels, resulting in a 39 × 39 correlation matrix. We zeroed all the negative correlation coefficients because of their ambiguous biological meanings [48,49]. We restricted following analyses to only positive correlations. All correlation coefficients were subjected to Fisher's z-transformation to improve normality.
Second, we binarized the correlation matrix according to sparsity-based threshold [50,51]. The correlation matrix was thresholded over a range of sparsity (from 20% to 40%, with 1% step size) in order to investigate the relationship between sparsity and the network properties. Specifically, each functional connectivity matrix Z ij can be converted to a binarized matrix B ij , where B ij is 1 if the value of the z value in matrix Z ij is greater than a given sparsity threshold and 0 otherwise [52]. Consistent with recent recommendations [52], we also adopted a single sparsity (25%) to normalize the network metrices, thereby allowing the comparison of group differences under the same topological organizations.

Network Analysis
In the current study, we used the GRETNA toolbox [53] in MATLAB 2022a to assess the topological measures, including global network metrics (small-world Gamma, Lambda, Sigma, global efficiency, and local efficiency) and regional nodal metrics (global nodal efficiency and local nodal efficiency) for the brain network.

Global Network Metrics
Small world. We focused on two key properties of the small world in a graph G: the cluster coefficient (C p ) and the characteristic path length (L p ) [54]. C p is the average of cluster coefficients over all nodes in the network. One-node cluster coefficient refers to the number of existing edges between that node and its neighbors, divided by the number of all theoretically possible edges [55]. For a given graph G with n nodes and K edges, the C p of the graph G is computed as follows [54]: where D i is the number of edges connected to node i. E i denotes the number of edges in the subgraph. C p reflects the local interconnectivity and cliquishness of a brain network [47].
The L p refers to the average of the shortest path lengths between all pairs of nodes in the network, whereas the shortest path length is defined as the minimum edges that link arbitrary nodes [55]. The L p of a graph G is defined as the average of the shortest path lengths between all pairs of nodes in network [52,54]: where d ij represents the shortest path length between node i and node j. Therefore, L p reflects the ability of serial information propagation within the network [55].
To examine the small-world attributes of a network, the normalized C p (referred to as Gamma) and normalized L p (referred to as Lambda) were calculated. Gamma and Lambda were computed as the ratio of the values to random rewired networks [56]: and L real p denote the clustering coefficient and the characteristic path length of a real network, whereas C rand p and L rand p are the means of the same parameters derived from 1000 matched random networks. The random network has the same number of nodes and edges and the same distribution of degrees as the real one. The ratio of Gamma to Lambda is defined as Sigma [57]: Gamma Lambda A small-world network is typically characterized as Lambda ≈ 1, Gamma > 1, and Sigma > 1 [58].
Efficiency. We also computed two properties of network regarding its efficiency, i.e., global efficiency and local efficiency. Global efficiency is defined as the inverse of the harmonic mean of the shortest path lengths between two arbitrary nodes in the entire network [55]. Global efficiency is calculated as follows: where d ij represents the shortest path length between node i and node j. Global efficiency reflects parallel information transformation at the global level in a network [50]. Local efficiency is computed as the average efficiencies of all nodes [55], wherein one given node and its direct neighbors compromised a sub-network. Local efficiency is calculated as follows: where E glob is the global efficiency of G i , which is the subgraph of the neighbors of node i. According to its definition, local efficiency can be regarded as a measure of information transfer within the immediate neighborhood of each node [50].

Regional Nodal Metrics
Complementary to global network metrics, we assessed nodal global efficiency and nodal local efficiency to provide measures for efficiency at regional level. For an indexed node, nodal global efficiency (also known as nodal efficiency) is defined as the inverse of the harmonic mean of minimum path length between that node and all other nodes in the network [50,59]: where d(i,j) denotes the shortest path length between node i and node j. It measures how well a sub-group is integrated in the whole network and reflects the ability of information exchange of the node itself [60].
The nodal local efficiency is measured as follows: where G i is the subgraph that includes node i and all its direct neighbors. Nodal local efficiency measures the communication ability of a sub-network consisting of the node itself and its direct neighbors [60].
The node of high efficiency is important for information integration and distribution [61]; in this study, the nodes with higher values in nodal efficiency (at least 1 SD larger than the average of all nodes in the brain network) were defined as brain hubs that are typically assumed to play critical roles in the functional integrity of whole networks [52]. BrainNet Viewer was used for visualization of regional nodal properties [62]. Channel-wise independent-sample t-tests were conducted to compare the group differences in nodal global efficiency and nodal local efficiency, with the false-discovery-rate (FDR) method accounting for the multiple comparisons [63].

Support Vector Regression
To explore whether it is possible to predict depression level (as assessed by HAMD scores) based on brain network metrices, we used the epsilon-support vector regression (ε-SVR, [64]). The radial basis function was utilized to construct the non-linear SVR model. We employed a grid search-based approach for hyperparameter optimization to determine the optimal regression parameters (i.e., C, γ, and ε). A nested cross-validation method was implemented [65], with the outer leave-one-out cross validation (LOOCV) estimating the generalization performance of the model and inner 10-fold CV estimating and selecting the optimal hyperparameters. The prediction accuracy was estimated using the Pearson correlation coefficient between the predicted and actual values [66]. The coefficient of determination (denoted by R 2 ) was also reported. We used the libsvm toolbox and custom codes in MATLAB to perform SVR analyses [67].

Resting-State Functional Connectivity
The RSFC pattern at the group level for patients with depression and healthy controls were illustrated in Figure 2. To explore whether there were significant differences in RSFC between the two groups, a series of independent-sample t-tests were conducted. Interestingly, the results showed that patients with depression (M ± SD, 0.32 ± 0.06, Pearson correlation coefficients) and healthy controls (0.31 ± 0.06) exhibited similar RSFC. Thus, RSFC analysis was insufficient to discriminate between the two groups in this study. correlation coefficients) and healthy controls (0.31 ± 0.06) exhibited similar RSFC. Thus, RSFC analysis was insufficient to discriminate between the two groups in this study.

Global Network Properties
To further investigate the intrinsic organization of cortical networks, we next analyzed global network properties, including small-world parameters (Gamma, Lambda, and Sigma) and global and local efficiencies, for the two groups. Figure 3 shows the profiles of five global network properties as functions of the sparsity thresholds (ranging from 20% to 40%, with 1% step size). We observed that the depression group showed lower small-world Lambda but higher global efficiency than the control group.

Global Network Properties
To further investigate the intrinsic organization of cortical networks, we next analyzed global network properties, including small-world parameters (Gamma, Lambda, and Sigma) and global and local efficiencies, for the two groups. Figure 3 shows the profiles of five global network properties as functions of the sparsity thresholds (ranging from 20% to 40%, with 1% step size). We observed that the depression group showed lower small-world Lambda but higher global efficiency than the control group.
Brain Sci. 2022, 12, x FOR PEER REVIEW 7 of 16 correlation coefficients) and healthy controls (0.31 ± 0.06) exhibited similar RSFC. Thus, RSFC analysis was insufficient to discriminate between the two groups in this study.

Global Network Properties
To further investigate the intrinsic organization of cortical networks, we next analyzed global network properties, including small-world parameters (Gamma, Lambda, and Sigma) and global and local efficiencies, for the two groups. Figure 3 shows the profiles of five global network properties as functions of the sparsity thresholds (ranging from 20% to 40%, with 1% step size). We observed that the depression group showed lower small-world Lambda but higher global efficiency than the control group. In accordance with recent studies [52], we adopted a single sparsity (i.e., 25%) to normalize all of the networks to explore the group differences in the same-size network topological organization. We measured lower Lambda in the depression group (1.12 ± 0.04), as relative to the control group (1.16 ± 0.10, t = 2.11, p = 0.04, Cohen's d = 0.53). We also found that the depression group (0.51 ± 0.03) showed higher global efficiency, compared to the control group (0.48 ± 0.05, t = 2.19, p = 0.03, Cohen's d = 0.73; Figure 4). In accordance with recent studies [52], we adopted a single sparsity (i.e., 25%) to normalize all of the networks to explore the group differences in the same-size network topological organization. We measured lower Lambda in the depression group (1.12 ± 0.04), as relative to the control group (1.16 ± 0.10, t = 2.11, p = 0.04, Cohen's d = 0.53). We also found that the depression group (0.51 ± 0.03) showed higher global efficiency, compared to the control group (0.48 ± 0.05, t = 2.19, p = 0.03, Cohen's d = 0.73; Figure 4).

Regional Nodal Properties
Having established the group differences in global network properties, we then tested the differences between the depression and control groups in terms of their regional nodal properties. Regarding nodal global efficiency, we found no hub in the depression group, but detected two central hubs (channels 26 and 30) in the control group ( Figure  5A,B). Concerning nodal local efficiency, we found three frontal hubs (channels 3, 15, and 16) and three central hubs (channels 25, 31, and 39) in the depression group and six central hubs (channels 22, 25, 28, 29, 30, and 33) and one occipital hub (channel 36) in the control group ( Figure 5C,D). Critically, the depression group (0.56 ± 0.36) compared to the control group (0.81 ± 0.20) showed a significantly lower nodal local efficiency at channel 36 (t = 3.38, pFDR < 0.05, Cohen's d = 0.86), which roughly corresponds to the left middle occipital gyrus [68]. As a robustness check, we additionally ran a one-way ANCOVA on the nodal local efficiency at channel 36, with HAMD scores and medication usage as covariates. The results showed that the group effect was still significant (F = 10.87, p = 0.0017, ηp 2 = 0.16). There were no significant effects for the covariates (F = 0.98, p = 0.37 for HAMD score and F = 0.24, p = 0.62 for medication usage). No significant group differences in nodal global efficiency were found (all pFDR > 0.05).

Regional Nodal Properties
Having established the group differences in global network properties, we then tested the differences between the depression and control groups in terms of their regional nodal properties. Regarding nodal global efficiency, we found no hub in the depression group, but detected two central hubs (channels 26 and 30) in the control group ( Figure 5A,B). Concerning nodal local efficiency, we found three frontal hubs (channels 3, 15, and 16) and three central hubs (channels 25, 31, and 39) in the depression group and six central hubs (channels 22, 25, 28, 29, 30, and 33) and one occipital hub (channel 36) in the control group ( Figure 5C,D). Critically, the depression group (0.56 ± 0.36) compared to the control group (0.81 ± 0.20) showed a significantly lower nodal local efficiency at channel 36 (t = 3.38, p FDR < 0.05, Cohen's d = 0.86), which roughly corresponds to the left middle occipital gyrus [68]. As a robustness check, we additionally ran a one-way ANCOVA on the nodal local efficiency at channel 36, with HAMD scores and medication usage as covariates. The results showed that the group effect was still significant (F = 10.87, p = 0.0017, η p 2 = 0.16). There were no significant effects for the covariates (F = 0.98, p = 0.37 for HAMD score and F = 0.24, p = 0.62 for medication usage). No significant group differences in nodal global efficiency were found (all p FDR > 0.05).

Prediction of the Depression Level
Finally, we sought to test whether it is possible to predict the depression level based on intrinsic organization of cortical networks in the depression group. To this end, we constructed ε-SVR models with nested cross-validation ( Figure 6A). We extracted global and local network properties, respectively, for each participant as predictors. The depression level, as measured by the HAMD scores, was the outcome variable. A data-driven method was used to search for the best predictor for the prediction analyses (one type of network properties as the predictor each time). We found that only when the occipital hub (i.e., channel 36) was derived from the nodal local efficiency and used as a predictor, was the correlation between actual and predicted HAMD scores significant (r = 0.45, p = 0.01, R 2 = 0.15) in the depression group ( Figure 6). These results indicate that it is possible to infer the depression level based on patients' nodal local efficiency in the occipital hub.

Prediction of the Depression Level
Finally, we sought to test whether it is possible to predict the depression level based on intrinsic organization of cortical networks in the depression group. To this end, we constructed ε-SVR models with nested cross-validation ( Figure 6A). We extracted global and local network properties, respectively, for each participant as predictors. The depression level, as measured by the HAMD scores, was the outcome variable. A data-driven method was used to search for the best predictor for the prediction analyses (one type of network properties as the predictor each time). We found that only when the occipital hub (i.e., channel 36) was derived from the nodal local efficiency and used as a predictor, was the correlation between actual and predicted HAMD scores significant (r = 0.45, p = 0.01, R 2 = 0.15) in the depression group ( Figure 6). These results indicate that it is possible to infer the depression level based on patients' nodal local efficiency in the occipital hub.

Discussion
In the present study, we performed functional connectivity-based comparisons between patients with depression and healthy controls. The results showed that patients with depression and healthy controls exhibited similar RSFC. However, the depression group showed lower small-world Lambda and higher global efficiency than the control group. The depression group compared to the control group showed a significantly lower nodal local efficiency at the occipital hub (i.e., channel 36), which roughly corresponds to

Discussion
In the present study, we performed functional connectivity-based comparisons between patients with depression and healthy controls. The results showed that patients with depression and healthy controls exhibited similar RSFC. However, the depression group showed lower small-world Lambda and higher global efficiency than the control group. The depression group compared to the control group showed a significantly lower nodal local efficiency at the occipital hub (i.e., channel 36), which roughly corresponds to the left middle occipital gyrus. When the occipital hub (channel 36) was derived from the nodal local efficiency and used as a predictor, the correlation between actual and predicted HAMD scores was significant in the depression group.
Resting-state functional MRI (rs-fMRI) revealed reduced nucleus accumbens functional connectivity in default mode network (DMN) in patients with recurrent major depressive disorder [69]. However, the decreased DMN functional connectivity was related with medication usage but not with MDD duration, as it was not found in first-episode drugnative MDD [70]. In this study, our results did not support our first hypothesis. The MDD patients were not drug-naïve, and it is possible that the drug usage in our patients might explain why their RSFC is similar to the healthy controls.
Our second hypothesis was verified. In this study, a more integrated brain network (lower Lambda and higher global efficiency) was found in MDD patients despite some negative results (Gamma, Sigma, and local efficiency). A review using non-invasive neuroimaging data and graph theoretical approaches for psychiatric disorders found that patients with depression did not display consistent alterations in small-world properties [71]. However, a significant increase in Gamma and Sigma despite the increase in Lambda was found in patients with MDD after electroconvulsive therapy [72]. It suggested that the Lambda may be a more stable trait biomarker of depression, despite the treatment. Some researchers detected a decreased global efficiency within MDD patients [39,73,74], which may be related to negative affective processing [75]. Other studies found no significant differences between MDD and healthy controls in terms of global efficiency [76,77], or generally higher global efficiency in MDD individuals [78]. Increased local efficiency, which associated with high HAMD score, was found in individuals with depression symptoms [79,80]. That indicated excessively high network segregation which might increase the brain's overall wiring cost [81]. Aforementioned studies were based on fMRI or EEG, while the present study was based on fNIRS, which found a different topological change in the brain network of MDD patients. It is likely that the metric of global efficiency might be sensitive to subtypes of depression, which deserves future investigations.
Crucially, this study found the MDD patients had significantly lower nodal local efficiency at the occipital gyrus that could predict HAMD scores. Our results echo a previous study that indicated aberrant nodal efficiency and centrality of regional connectivity in the occipital cortex [82]. Regional homogeneity (ReHo) was found to be lower in the occipital gyrus among MDD patients, which could even help to discriminate patients with melancholic MDD from patients with non-melancholic MDD [83]. In addition, MDD patients had abnormal local intrinsic gray-matter connectivity in the occipital cortex, which was associated with some symptoms of depression [84]. These results might help explain the aberrant topological properties of brain functional connectivity at the occipital hub we found in MDD patients. Researchers found a smaller gray matter volume in the left occipital middle gyrus in MDD patients, compared with controls [85,86]. As for white matter, fractional anisotropy in the left middle occipital gyrus was reduced in MDD patients [87]. Meanwhile, MDD patients had decreased cerebral blood flow in the left middle occipital gyrus, especially in those with acute phase and medication-free [88]. These results might be a structural basis of the lower nodal local efficiency at the left middle occipital gyrus. However, the abnormal activity in the left middle occipital gyrus may be state-specific in current and remitted MDD patients [89]. That may help explain some inconsistent results [90].
Our findings at the occipital hub could also be interpreted from an "information processing" perspective [91]. A recent study found that the amplitude of low-frequency fluctuations in the left middle occipital gyrus decreased within patients with MDD, compared to healthy controls; the researchers argued that the left middle occipital gyrus may be involved in the processing of cognitive biases of MDD in resting states [92]. Consistent with this research, decreased nodal local efficiency at the middle occipital gyrus bolstered the notion that processing bias in MDD may be initiated as a perceptual visual bias. Biased information from the occipital hub could then spread through the brain due to the higher global efficiency in MDD patients. As the DMN hyperactivity was related to negative rumination in depression [93], the biased information initiated from the occipital hub may, eventually, cause a series of cognitive and affective symptoms of MDD.

Limitations and Strengths
There were some limitations in our study. Firstly, the MDD patients were not medicationnaïve in this study. Thus, we could not exclude the possibility that these drugs affected the current results. Secondly, the sample was inadequate in looking for a differential effect across depressive subtypes. Thirdly, fNIRS has a penetration limit into the superficial gray matter of the cortex of around 2 cm. Therefore, we could not detect subcortical changes. Nevertheless, the topology of the cortex in MDD patients did show some abnormities. The main strength of this study is that we used a more convenient way of rs-fNIRS combined with graph theory to distinguish patients with MDD from healthy people.

Conclusions
In conclusion, we found a more integrated brain network with lower Lambda and higher global efficiency in MDD patients. Meanwhile, MDD patients had a lower nodal local efficiency at the occipital hub, which could predict depressive symptoms. These results may provide a new method in assisting the clinical diagnosis and help to elucidate the brain mechanism of MDD. Future studies could verify these results in subgroups of MDD patients, such as untreated patients, patients in remission stage, male and female patients, young and old patients, and so on. In addition, the other regions or areas should be explored.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to restrictions of privacy and ethic.