Frequency Dependant Topological Alterations of Intrinsic Functional Connectome in Major Depressive Disorder

Major depressive disorder is associated with aberrant topological organizations of brain networks. However, whether this aberrance is shown in broader frequency bands or in a specific frequency band remains unknown. Fifty patients and fifty gender, age and education matched normal controls underwent resting state functional magnetic resonance imaging. Frequency dependent topological measures based on graph theory were calculated from wavelet decomposed resting state functional brain signals. In the specific frequency band of 0.03–0.06 14Hz, the clustering coefficient and the global efficiency were reduced while the characteristic path length was increased. Furthermore, patients showed aberrant nodal centralities in the default mode network, executive network and occipital network. Network based statistical analysis revealed system-wise topological alterations in these networks. The finding provides the first systematic evidence that depression is associated with frequency specific global and local topological disruptions and highlights the importance of frequency information in investigating major depressive disorders.

and amnesic mild cognitive impairment 25 . With respect to MDD, only one study used wavelet analysis to investigate brain network configurations in MDD patients 28 . However, this study by Manoliu and his colleagues 28 only investigated a specific frequency band of 0.060 , 0.125 Hz, which may have overlooked valuable information in other frequency bands.
Here we measured functional connectivity and examined topological organization of MDD patients in resting state fMRI (R-fMRI) data. In brief, processing procedures included (1) constructing network matrices based on wavelet decomposed R-fMRI data, (2) calculating network topological metrics, (3) comparing topological metrics across groups for each wavelet scale and (4) correlating topological metrics with clinical variables. Given previous graph theory based studies of MDD [17][18][19][20]28 , We hypothesized that the brain networks of MDD patients would be disrupted at both global and regional levels. A recent study showed that intelligence was mediated by the coupling of brain networks across several frequency bands 29 . We thus also hypothesized that the aberrant brain networks of MDD patients may be manifestated across multiple frequency intervals.

Results
Demographic statistics. The depression group and healthy control (HC) group were matched in gender, age and education (all p . 0.05). The two groups significantly differed in HRSD  Table 1.
Frequency specific global topological alterations in patients. Significant group differences in global topology were exclusively found in Scale 3 (0.03-0.06 Hz). First, in the specific thresholded network of Scale 3, patients and HCs differed significantly in the total number of connections and the mean correlation as well as the mean anatomical distance (mean Euclidean distance across existing edges). The patient group featured less connections, weaker connections and shorter anatomical connections (mathematically defined as the Euclidean distance between centroids of every two nodes) in Scale 3 (See Figure 1). Second and further, as shown in Figure 2, global topological metrics exhibited significant alterations in Scale 3 but not in the other Scales. In patients, the clustering coefficient and the global efficiency were reduced while the characteristic path length was increased.
Aberrant nodal centrality. Degree centrality was altered (p , 0.05 uncorrected) in several cortical and subcortical regions (see Figure 3 and Table 2). Patients had decreased degree centrality (i.e. less connections) predominately in nodes relevant to executive function, emotion processing and basic sensory function. They also showed increased betweenness centrality in cortical midline structures (CMS), which are part of the default mode network and are a core component of self processing.
Aberrant functional connectivity. For a liberal primary threshold of p 5 0.005, the NBS method detected a dysconnected network (p , 0.05. corrected) comprised of 136 decreased connections linking wide spread regions (see Figure 4A, Table 3 lists the nodes). Under a strict threshold of p 5 0.001, the NBS method indentified a smaller network comprised of 23 reduced functional connections in patients (see Figure 4B, Table 3 lists the nodes). Of note, networks identified by the NBS method can only be statistically significant as a whole. Any connections within cannot be declared as independently significant 30 .

Correlation between topological metrics and clinical variables.
None of the global topological metrics significantly covaried with HRSD score or disease duration. Among the regions showing significant group differences, the bilateral insular cortex and  bilateral planum polare exhibited trend of positive correlations with HRSD scores (see Figure 5), but they did not survive multiple corrections.

Discussion
The present study examined the topological organization of the brain networks of MDD patients in R-fMRI data and found topological alterations in the specific and most salient wavelet scale of 0.03 , 0.06 Hz 14 . Globally, we observed a reduction in the clustering coefficient and global efficiency as well as an increase in characteristic path lengths. Regionally, we found aberrant nodal centralities in widely distributed cortical and subcortical regions. Finally, the NBS method identified dysconnected subsystems (under primary thresholds of 0.005 and 0.001 respectively) in patients. Significant global topological differences between networks of MDD patients and HCs were only found in the specific frequency band of 0.03 , 0.06 Hz. Previous researchers have suggested that neuronal oscillations are linearly arranged on the natural logarithmic scale and this regularity indicates that different frequency bands result from different oscillators with different biological properties and functions 21,22 . Although the exact mechanisms underlying signals of different frequency bands remains poorly understood, the 0.03 , 0.06 Hz range proved to be most salient resting brain state in normal adults 14 . Moreover, disrupted topological architectures of amnestic mild cognitive impairment (aMCI) patients 25 and Alzheimer's disease (AD) patients 24 as well as regional alterations of Attention deficit hyperactivity disorder (ADHD) patients 23 were in this specific frequency band. A recent work highlight intelligence of healthy adults was mediated by coupling between multiple networks across multiple frequencies 29 , This evidence in conjunction with the present results and previous studies might serve to highlight the clinical significance of brain oscillations under 0.03 , 0.06 Hz.
The human brain functions like a small world 14 . In the powerful framework of graph theory, a network is declared as small world when it attains a perfect balance between local specialization (indexed by a high clustering coefficient) and global integration (indexed by low characteristic path length) 13 . The present study found the brain networks of MDD patients and HCs both exhibited the organization principles of small worldness.
Though networks of MDD patients preserve small worldness, the configurations of networks in MDD patients were perturbed. First, networks of MDD patients demonstrated fewer connections and weaker connections. Compatible with our finding in the R-fMRI data, studies based on DTI data also indicated sparse and reduced white matter connectivity in MDD patients 31 . Second, compared to HCs, the clustering coefficient and global efficiency of networks in MDD patients were reduced while characteristic path length was increased. This suggests that the architecture of networks in patients were less locally specialized and less globally integrated.
Degree centrality was decreased in 38 regions of MMD patients' brain networks. Most of these regions were relevant to cognitive executive functions, emotion processing and basic sensory functions. Within the executive control network (which mainly involves the lateral frontal and superior parietal lobules) 32 , the present study indentified fewer global connections in the frontal pole, inferior frontal gyrus and posterior superamarginal gyrus. This may suggest that patients were less skilled in executive functioning as supported by working memory 33 and Go/Nogo studies 34 . With respect to emotion processing, researchers have documented impaired prefrontallimbic circuits in patients with MDD 35 . Compatible with these results, we found decreased degree centrality in the prefrontal cortex and limbic areas. Specifically, some researchers have underlined the pathophysiological role of the subcallosal cingulate gyrus and amygdala in depression 36,37 . The observed reduced degree centrality in the prefrontal-limbic circuit may be attributable to inadequate reciprocal interactions between the prefrontal cortex and limbic areas. Finally, patients with MDD were found to have low concentrations of caminobutyric acid (GABA) 38 and efficient treatment brings GABA to presymptomatic levels 39 . The observed reduction of occipital regions may result from disruptions to the metabolism of the occipital cortex in MDD patients.
Betweenness centrality was increased mainly in cortical midline structures (CMS) 40 . CMS refer to brain structures situated near the medial wall of the brain and is mainly comprised of the ventral medial prefrontal cortex, dorsal medial prefrontal cortex, anterior cingulate cortex (the supragenual part in particular) and posterior cingulate cortex. CMS are core biological components for the self 40 . Depression is characterized by a stronger tendency toward introspecting and self-focused cognitions 41 . Hyperactivity in CMS has been recorded in MDD patients 42 . Therefore, the elevated betweenness centrality in patients may reflect the strengthened communication among CMS regions and constitute the pathophysiological basis of rumination in MDD patients.
The NBS method was implemented to locate aberrant connections in MDD patients. A larger disconnected network was identified under the liberal threshold (p 5 0.005). Most nodes of this network were those nodes with altered nodal centrality. Furthermore, when we employed a strict threshold (p 5 0.001), a smaller abnormal network featuring a majority of long-range connections was found. This suggests that deficient global information integration in patients was mainly driven by decreased long-range connectivity. Of note, this smaller network included nodes within the prefrontal-limbic circuit and provides empirical support for the proposed pivotal role of the prefrontal-limbic circuit in the pathophysiology of depression [35][36][37] .
Of additional note, the disease duration of the patients here were above 4 years, which means that the current finding may reflect the chronic aspects of depression disorder. Previous psychological studies already indicated that chronically depressed patients differed from acute depressed patients in cognitive functions 43 , response to medications 44 . However, we find few studies investigated the aberrant brain mechanisms of medication free chronic depression. Therefore, the present finding may help gain insight into this particular field and future neuroimaging studies which compared chronically depressed patients with acute ones may further shed lights to the specific brain dysfunction of chronically ones.
The present study is not without limits. First, patients in our study included 17 medicated patients. Some researchers reported a normalization influence of depression treatments on functional connectivity 45,46 . However, the exact mechanism remains largely unknown 47 . Moreover, we think the present results were not likely to have been driven by medication effects for the following reasons: (1) medication status (0 or 1) did not significantly covary with any global topological metrics (all p . 0.05); (2) global topological metrics did not significantly differ between medicated and nonmedicated patients; (3) when we regressed out medication status, similar results were obtained (i.e. significant global topological differences were exclusively found in wavelet scale. Second, the significances of group comparisons in nodal centrality did not survive stringent bonferroni or liberal fdr corrections. Therefore, the results should be viewed as an explorative attempt and could be reaffirmed by methods with more statistical power (larger samples or selecting prior regions of interest). Third, methodological selections such as parcellation schemes, threshoding methods, smoothing effects, network types (binary or weighted) and modality types may have impacted upon our findings 48 . This may partly account for the discrepancy between finding across network based studies of MDD patients [17][18][19][20]28 .
To our knowledge, our work is the first to explore frequency dependent topological alterations in major depressive disorder. Compared to age, gender and education level matched controls, patients showed both global and regional topological alterations in the specific and most salient frequency range of 0.03 , 0.06 Hz. The finding provides the first systematic evidence that depression is associated with frequency specific global and local topological disruptions and highlights the importance of frequency information in investigating major depressive disorders.

Methods
Participants. Fifty patients diagnosed with major depressive disorder and fifty normal controls participated in this investigation. No participants were left-handed as measured by the Edinburgh Handedness Inventory 49 . Written informed consent  50 . All patients had a minimun HRSD score of 18. Nineteen patients had been receiving antidepressant therapy, but were free of any treatment 4 weeks before recruiment. Normal controls were recuited from the local community around Southwest University, Chongqing. All controls scored less than 7 on the HRSD. No participants met the exlusion criteria of any neurological or mental disorder besides depreesion, substance abuse or head trauma. The principls of the present study was adequately guided the approved guidelines which was in accordance with the Decleartion of Helsinki. The study was approved by the Southwest University Brain Imaging Center Institutional Review Board.
Image Acquisition. All images were acquired on a 3.0-T Siemens Trio MRI scanner using a 32-channel whole-brain coil (Siemens Medical, Erlangen, Germany). Highresolution T1-weighted 3D images were acquired using a magnetization-prepared rapid gradient echo (MPRAGE) sequence (echo time . During image acquisition, participants were instructed to keep their eyes closed while keeping their head as still as possible without falling asleep. All participants stayed awake during MRI examination which was confirmed by the participant after the examination. Data preprocessing. Data preprocessing was performed using DPARSF (http:// restfmri.net). The first 10 EPI scans were discarded to suppress equilibration effects. The remaining 232 scans of each subject underwent slice timing correction by sinc interpolating volume slices, motion correction for volume to volume displacement, spatial normalization to standard Montreal Neurological Institute (MNI) space using affine transformation and nonlinear deformation with a voxel size of 3 3 3 3 3 mm 3 , temporal high pass filtering ($0.01 Hz) to reduce the effect of low frequency drift, and regressing out 24 head-motion parameters (6 motion parameters for current the volume, 6 motion parameters for the previous volume and 12 corresponding squared items). Originally proposed by Friston and his colleagues 51 , the Friston 24 parameter model outperformed other motion correction strategies in reducing motion related effects in a recent work 52 . Furthermore, no subjects were found to have excess head motion (translation . 1 mm or rotation . 1) and head motion profiles were statistically matched in any direction (all p . 0.05) between groups. Given the debate over removing global signal in R-fMRI data preprocessing 53,54 , we did not regress the global signal as had been done in previous wavelet studies 14,24,27 . To exclude cofounding effect from CSF and white matter signals, we further regressed out the CSF and WM signals. A custom grey matter mask was created by thresholding the mean GM probability map from all individuals (threshold 5 0.2). This procedure was implemented in SPM 8 (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/). All subsequent analyses were restricted to the custom grey matter mask.
Network construction. Node definition. Nodes were obtained by anatomically dividing the brain into 112 distinct regions according to the Harvard-Oxford atlas 55 .
Edge definition. First, the mean time series was extracted from each region. Second, the maximal overlap discrete wavelet transform (MODWT) method 56 was used to decompose each regional mean R-fMRI time series into four successive scales or frequency intervals ( 59 . Second, for the selected wavelet scales, nodal centrality was calculated for each node and compared across groups. Here, nodal centrality refers to degree centrality and betweenness centrality 59 . Regional analyses were exploratory to get a complete understanding of nodal aberrancies. All topological metrics were automatically computed by GRETNA (https://www. nitrc.org/projects/gretna/). Detailed definitions of topological metrics can be found in the Supplementary Materials section.
Statistical analysis. Group comparisons. Group comparisons of topological metrics (global and regional) were carried out by permutation tests (10000 permutations) while controlling for gender, age and education. Briefly, for each permutation, subjects were randomly shuffled and the group difference of a given metric was recomputed. We obtained a null distribution of group differences of that metric and the significance level of the real group difference was determined by the percentage of permutations with a group difference equal to or greater than the real group difference in all permutations (10000). In order to find aberrant connections in patients, we carried out a method called network based statistical (NBS) analysis 30 . Specifically, first, edges with a p value (obtained by a two sample t-test) less than a primary uncorrected threshold (p , 0.005 and p , 0.001 respectively) were identified. Based on these superathreshold edges, several components or subnetworks were then identified and their sizes (number of edges within a subnetwork) calculated. Second, the significance level of each identified components was computed according to a null distribution of the largest subnetwork size which was empirically obtained through 10000 permutations. The significance level of any originally found subnetwork (i.e. before permutation test) with size M was computed as the percentage of permutations in which the size of the largest component exceeded or equaled M. Only significant subnetworks (p , 0.05, family wise error corrected) were reported. The reason we use two primary threshold (here, p , 0.005 and p , 0.001 respectively) is to give a full understanding of aberrant connections in patients as done in previous studies 25, 60 . In addition, we use fdr (p , 0.05 corrected) here for multiple comparisons (2 comparisons here: p , 0.005 and p , 0.001).
Correlation analysis between topological metrics and clinical variables. Partial correlation analysis was employed to assess the relationships between topological metrics and clinical variables (HRSD scores and disease durations) while regressing out nuisance variables (i.e. gender, age, education).