Effect of SSRIs on Resting-State Functional Brain Networks in Adolescents with Major Depressive Disorder

Investigation of brain changes in functional connectivity and functional network topology from receiving 8-week selective serotonin reuptake inhibitor (SSRI) treatments is conducted in 12 unmedicated adolescents with major depressive disorder (MDD) by using wavelet-filtered resting-state functional magnetic resonance imaging (fMRI). Changes are observed in frontal-limbic, temporal, and default mode networks. In particular, topological analysis shows, at the global scale and in the 0.12–0.25 Hz band, that the normalized clustering coefficient and smallworldness of brain networks decreased after treatment. Regional changes in clustering coefficient and efficiency were observed in the bilateral caudal middle frontal gyrus, rostral middle frontal gyrus, superior temporal gyrus, left pars triangularis, putamen, and right superior frontal gyrus. Furthermore, changes of nodal centrality and changes of connectivity associated with these frontal and temporal regions confirm the global topological alternations. Moreover, frequency dependence is observed from FDR-controlled subnetworks for the limbic-cortical connectivity change. In the high-frequency band, the altered connections involve mostly frontal regions, while the altered connections in the low-frequency bands spread to parietal and temporal areas. Due to the limitation of small sample sizes and lack of placebo control, these preliminary findings require confirmation with future work using larger samples. Confirmation of biomarkers associated with treatment could suggest potential avenues for clinical applications such as tracking treatment response and neurobiologically informed treatment optimization.


Introduction
Major depressive disorder (MDD) is a highly debilitating condition that can result in tragic outcomes such as chronic disability and suicide. The onset of depression is frequent during the adolescent period [1], a time notable for significant brain development and physical and social changes [2,3]. According to the World Health Organization (WHO), more than 10% of adolescents are affected by MDD in the US [4,5]. Moreover, early onset MDD has been shown to increase the risk of developing adult depression [6]. However, current evidence-based treatments including antidepressant medication and cognitive behavioral therapy are only successful in reducing depression in about half to two-thirds of cases [7]. To design and optimize new and more effective treatments, a deeper understanding of the neurobiological basis of existing treatments is needed.
During the last decade, MRI has been extensively applied in studying neural underpinnings of depression. Using structural MRI, previous studies have reported that patients with depression show reduced volumes compared to healthy controls in the anterior cingulate cortex, orbitofrontal cortex, and hippocampus [8][9][10]. Based on diffusion MRI, previous studies have shown disrupted white matter integrity in the limbic system, dorsolateral prefrontal cortex, thalamic fibers, and corpus callosum [11][12][13][14][15]. Additionally, numerous reports have shown disruptions of functional connectivity in the limbic system using functional MRI (fMRI) [16][17][18][19][20]. Taken together, these findings suggest that MDD involves a complex set of connectivity deficits in the fronto-limbic system as well as in other brain networks.
Since MDD involves neural network abnormalities, a growing body of research has been dedicated to understanding how such networks may change after antidepressant treatments [21][22][23][24][25][26]. However, prior studies have largely focused on specific regions of interest, i.e., regional functional activation and functional connections from them, such as the amygdala. For example, in a prior report from the same dataset, as reported here, we investigated treatment-related changes in resting-state functional connectivity between the amygdala and all other voxels in the brain [21,22]. However, it is possible to examine connectivity changes at a more global level in order to assess change in connections among all brain regions. Furthermore, beyond functional connectivity analyses, a more advanced approach to understanding brain networks is to characterize the topological characteristics of the network. Brain network studies have recently shown the promising power of graph-theory based analysis in examining topological properties for brain networks considering cortical/subcortical regions as nodes and anatomical or functional connectivity among regions as edges [27][28][29][30][31]. Using graphical analysis, the small-world property has been discovered in both structural and functional networks for healthy brains [32][33][34][35][36]. A small-world network has the property of a high clustering coefficient and a low characteristic path length. These features are associated with efficient local information delivery, i.e., functional segregation, and effective collaboration for distributed information processing, i.e., functional integration, respectively, of the brain network structure [27]. Graph theory analysis has been extensively used to uncover brain network abnormalities in patients with schizophrenia [37][38][39][40], Alzheimer's disease [41][42][43][44], and depression [14,[45][46][47]. These promising results motivate the exploration of treatment-related topological patterns in this work from the rapidly changing adolescent brain networks, which has not been investigated in previous studies. Given the complexity of MDD symptoms and response to treatment, knowledge about possible topographical disruptions in functional brain networks could provide greater direction in treatment selection and development.
This study aims to discover treatment-related changes in brain functional network topology and whole-brain connectivity for 12 unmedicated adolescents with MDD before and after receiving an 8-week selective serotonin reuptake inhibitor (SSRI) medication treatment. We specifically analyzed brain networks with frequency-selective restingstate functional connectivity (RSFC) by using filtered fMRI signals in four frequency bands. First, we examine treatment-related changes in a set of topological measurements. Second, we examine treatment-related changes in functional connectivity between pairs of all brain regions. Third, we examine treatment-related change in subnetworks that are composed of the significantly altered connections. Finally, we examine correlations between clinical assessment and brain network features (topological measures and connectivity). The framework for data analysis is shown in Figure 1.

Participants
Participants of this study consisting of twelve adolescents aged 12 to 19 years diagnosed primarily with MDD were unmedicated before the baseline scan. They were invited to return for a second neuroimaging scan after they completed 8 weeks of antidepressant treatment under the care of their own provider. They represent a subset of a larger, cross-sectional neuroimaging study [19]. After the informed consent process, all participants completed a comprehensive diagnostic assessment. Interviews were conducted separately with adolescents and parents, and included Kiddie Schedule for Affective Disorders and Schizophrenia-Present and Lifetime Version [48] and the Children's Depression Rating Scales-Revised (CDRS-R) [49]. Self-report measures assessing symptoms in the past two weeks included the Beck Depression Inventory II (BDI-II) [50,51] and the Inventory of Depression and Anxiety Symptoms (IDAS) [52][53][54]. The IDAS provides a score for the following symptom dimensions: general depression, dysphoria, lassitude, insomnia, suicidality, appetite loss, appetite gain, ill temper, well-being, social anxiety, panic, and traumatic intrusion. Patients with any of the following exclusionary conditions were excluded: current use of psychotropic medications, intellectual disability, pervasive developmental disorder, substance use disorder, bipolar disorder, and schizophrenia. The demographic information is summarized in Table 1. The experimental procedures involving human subjects described in this paper were approved by the University of Minnesota Institutional Review Board in 2008 (case number 0804S30542).  [50]; * from the onset of the first episode, which may not be the current episode.

Pre-and Post-Treatment Clinical Assessment
Diagnosis is based on the Kiddie Schedule for Affective Disorders and Schizophrenia-Present and Lifetime version (K-SADS-PL) [48] interviews completed separately with adolescents and parents [21,22,55] and the Beck Depression Inventory-II (BDI-II) [50]. The BDI-II was also administered on the day of the post-treatment scan. Change in depression symptoms was calculated by subtracting the pre-treatment scores from the post-treatment BDI-II scores (negative scores represent improvement as indexed by a drop in depressive symptoms).

Image Acquisition and Processing
MRI scans are conducted by using a Siemens 3 Tesla TIM Trio scanner located at the Center for Magnetic Resonance Research, University of Minnesota. A T1-weighted highresolution magnetization prepared gradient echo sequence scan is obtained with following parameters: repetition time = 2530 ms, echo time = 3.65 ms, inversion time = 1100 ms, flip angle = 7 degrees, 1 mm slices, field of view = 256 mm, voxel size 1 × 1 × 1 mm 3 , and GRAPPA = 2.5 min. For the resting-state fMRI (rsfMRI), functional data were acquired by using an echo planar imaging sequence (EPI) with the following: 180 T2-weighted wholebrain functional volumes with 34 interleaved slices of 4 mm thickness and no skip, tilted to AC-PC alignment, repetition time = 2000 ms, echo time = 30 ms, flip angle = 90 degrees, field of view = 220 mm, and 64 × 64 matrix (voxel size 3.4375 × 3.4375 × 4 mm 3 ). Participants are asked to close their eyes and stay awake. Physiological data (respiration and cardiac traces) are collected during the entire scan.
Brain parcellation is performed using Freesurfer (http://surfer.nmr.mgh.harvard. edu, accessed 20 August 2018) following the Desikan-Killiany atlas to create 76 cortical/subcortical regions listed in Table 2 [56]. FreeSurfer output was visually inspected; identified errors were manually corrected on a slice-by-slice basis. The parcellation is registered to the fMRI data using FLIRT and FNIRT tools in FSL (Functional Magnetic Resonance Imaging of the Brain Software Library; http://www.fmrib.ox.ac.uk/fsl, accessed 20 August 2018) [57][58][59]. Functional images are preprocessed for motion correction, brain extraction, RETROICOR processing (regression of physiological signal, i.e., heart beat and respiration, which were recorded during the scan [60]), geometric distortion correction using the field map, and removing regression of signal from CSF and white matter. Data scrubbing was performed following the method of Power and colleagues [61], excluding any volume with a value for the temporal derivative of time courses' root mean squared head motion variance exceeding 8 and/or a framewise-dependent value exceeding 0.5, along with the previous volume and the 2 following volumes.

Resting-State Functional Connectivity and Network Construction
Previous studies have found that the functional connectivity between brain regions is frequency dependent [32], and disease-related alterations in functional connectivity are prone to specific frequency bands [44,62,63]. Therefore, in this study, the fMRI signal is filtered using a 4-level stationary discrete wavelet transform (SDWT) [64,65] with 'db4' wavelet into different frequency bands. Unlike traditional decimated wavelet transform, the SDWT satisfies translation-invariance. In SDWT, the downsamplers and upsamplers are removed, and the filter coefficients are upsampled. The wavelet scales approximately correspond to frequency ranges of 0.12-0.25 Hz, 0.06-0.12 Hz, 0.03-0.06 Hz, and 0.015-0.03 Hz, respectively. After SDWT filtering, functional connectivity between brain regions is computed based on the Pearson correlation between the average fMRI time courses at each frequency band separately for each subject. A frequency-specific 76-by-76 undirected graph is constructed based on the 2850 correlation coefficients for each frequency band.

Graph-Theoretic Analysis
Adaptive thresholding is employed to create a series of brain networks corresponding to graph densities from 10% to 50% for network topology analysis. Individual bias is eliminated, and the dynamic range of connectivity is normalized by adaptive thresholding. The core network has the smallest density and contains the top 10% strongest connections; by contrast, the most robust network has the largest density and contains the top 50% strongest connections. The range of densities is chosen to ensure the capability of the small-world measure [66].

Graph Measures
Graph structures are quantified by using graph-based measures including clustering coefficient, characteristic path length, smallworldness, global/local efficiency, participation coefficient, within-module-degree z-score, degree, strength, and betweenness centrality. A complete review of topological brain network measures can be found in the literature [29,67]. Symbols are defined as the following:

Clustering Coefficient
The clustering coefficient, C i , is a measure of the degree to which nodes in a neighborhood tend to cluster together, and it is defined as the ratio of the number of closed triangles to the number of connected triples of vertices in the two-hop neighborhood of a node, i. The higher the clustering coefficient is, the more tightly the node is connected to its neighbors. Therefore, a high local clustering coefficient potentially implies high local efficiency.

Characteristic Path Length
The characteristic path length L net quantifies the integration ability of a network. The definition of L net is the average distance between any two nodes in the network. To avoid the disconnection problem, i.e., the distance between some nodes is infinity, the harmonic mean version of the original definition is used in this study: where n is the number of nodes in N and d ij is the shortest path length between node i and node j.

Smallworldness
The small-world properties of a network are quantified by the clustering coefficient and the characteristic path length of the network. The clustering coefficient of the entire network C net is the mean C i of all nodes in the network. This global measure quantifies the cliquishness of a network. In order to diagnose the small-world properties, global clustering coefficients and characteristic path lengths are normalized by the same metrics estimated from random networks with the same number of nodes, edges, and degree distribution. The normalized network clustering coefficient is defined as C norm = C net /C rand , the normalized characteristic path length is defined as L norm = L net /L rand , and the smallworldness of the network is defined as C norm /L norm . A small-world network is expected to have high local clustering and low mean path length; therefore, the smallworldness is greater than 1 for a small-world network.

Degree
Degree, defined as k i = ∑ j∈N a ij , measures how well a node is connected to the other nodes, which can also be interpreted as the (nodal) centrality of the node. The degree is defined as the number of binary connections to the node. Therefore, a node with a large nodal degree might be considered a hub.

Local Efficiency
Local efficiency, E i , is a measure of how efficiently nodes exchange information in the one-hop neighborhood of node i, and it is positively related to the clustering coefficient.

Modules
A module, M, also known as a cluster, is defined as a non-overlapping division of a network that represents the structural segregation in the brain. In a cluster, nodes are abundantly connected, and nodes in different clusters are connected through a few hub regions (nodes). The interpretation of the module formation is that different structural modules are segregated by biological and functional characteristics.

Participation Coefficient
Participation coefficient, where k i (m) = ∑ j∈m a ij , assesses the diversity of inter-modular interconnections of individual nodes. The parameter k i (m) denotes the within-module degree of node i. In the brain network, a provincial hub, which is an important part in the facilitation of modular segregation, will have a high within-module degree but a low participation coefficient. Furthermore, nodes having both high within-module degree and participation coefficient facilitate global inter-modular integration as connector hubs.

Betweenness Centrality
The betweenness centrality is a measure of centrality defined by the fraction of all shortest paths in the network that contain a given node. Nodes with high values of betweenness centrality participate in a large number of shortest paths.

Graph Measure AUCs
For binary graph measures, we considered the area under curve (AUC) of the measure across a range of graph densities. When the range of graph densities is uniformly sampled, the AUC is just the average measure across graph densities.

Statistical Analyses
In order to assess the significance of a change, D, in brain network measures from baseline to post-treatment, we first estimated the empirical distribution of the null hypothesis, H 0 : no significant change due to the treatment, by performing a random sign permutation test. During the test, the feature differences between before and after treatment are randomly assigned with a positive or negative sign to calculate the average difference, D , for a realization of the random process. To test the hypothesis, 1,000,000 realizations are implemented. The probability that the average difference, D , is greater or smaller (depending on the sign of D) than the original measured average difference, D, is the empirical tail probability, i.e., p-value. With N realizations, the best achievable resolution of the p-value is 1 N , which is 10 −6 here. In other words, p-value differences smaller then 10 −6 are not distinguishable, i.e., estimation may be zero for any p-value less than 10 −6 . Although, the best resolution for random sign permutation tests is set to 10 −6 , the resolution is still limited by the sample size. In this study, 12 subjects yield 2 12 combinations and result in the best p-value resolution of 2 −12 = 0.0002441.
As an exploratory study, uncorrected p-values are reported, and minimal multiplecomparison correction (corrected for direct comparisons, i.e., by number of regions for regionrelated metrics and by number of connections for connection-related metrics) is performed with Bonferroni correction to highlight potential changes for 5% false-discovery rate.
In addition to the p-value, effect size and power are also reported. The effect size of the treatment-related change was measured using Cohen's d with pooled standard deviation [68]. The statistical power of the test with significance level 0.05 is calculated based on the mean of changes, the standard deviation of changes, and sample size (N = 12), and it is tested against the null hypothesis (µ 0 = 0) [69,70].
In order to assess possible correlations between neural network changes and clinical changes, we conduct regression analysis using 'fitlm' (Matlab) by taking the changes in brain network measures as main regressors and change in total BDI-II scores as the outcome variable. Additionally, to discover whether baseline network measures might serve as predictors of treatment response, we also perform regression analysis using beforetreatment brain network measures as main regressors, again with the change in total BDI-II scores as the outcome variable.
Previous research suggests that analyzing a subnetwork can effectively reduce the number of hypothesis tests for controlling the false positive rate and increase statistical power in connectomic analysis [71,72]. Here, we create a subnetwork consisting of significantly changed connections selected with a false discovery rate (FDR) at level α = 5% using Benjamini-Hochberg procedure [73] for each frequency band. In order to examine the broadness of treatment effect to a region in communicating with the other regions, withinsubnetwork nodal degree is calculated, which is defined by the number of connections connected to the node in the subnetwork. Figure 2 shows normalized clustering coefficient and smallworldness in mean curves across patients as functions of network densities for baseline and after treatment in all four frequency bands. Curves show results similar to previous studies in the small-world topology that the normalized clustering coefficient is greater than one and approaches one when the network density increases and the smallworldness decreases when the network density increases [32,33,35,38,42]. However, in addition to previous findings, results indicate the clustering coefficient is increased after treatment in all frequency bands for small network densities, i.e., core networks. Given no change in characteristic path length, smallworldness shows consistent change with the clustering coefficient. Regular networks have large clustering coefficient and low smallworldness, and random networks have low clustering coefficient and low smallworldness. The increased clustering coefficient and smallworldness indicate that core network topology has changed toward a small world network that is commonly observed from human brain networks. In terms of statistical significance, the AUC of the normalized clustering coefficient differences has p-value ranging from 0.083 to 0.47, and the smallest value is achieved in the 0.12-0.25 Hz band with an effect size of 0.57. For the smallworldness, the p-value ranges from 0.143 to 0.431, and the smallest value is achieved in the 0.12-0.25 Hz band with an effect size of 0.4765.

Local Topological Metrics
Next, the local network measures associated with each brain region, i.e., network nodes, are analyzed. The results show that treatment-induced topology changes are frequency-specific. After treatment, MDD patients show (p-value < 0.00066 = 0.05/76, Bonferroni correction for 76 regions) decreased local efficiency and clustering coefficient in right caudal middle frontal gyrus in the 0.12-0.25 Hz band. In the 0.06-0.12 Hz band, decreased participation coefficient is observed in left pars triangularis. In the 0.03-0.06 Hz band, there are decreases in betweenness centrality in the right precentral gyrus and in within-module-degree z-score in the left lateral occipital cortex. The corresponding pvalues, effect sizes, and power are listed in Table 3, and the boxplots of the changes and measurements before and after treatment are shown in Figure 3. No significant differences were found in the 0.015-0.03 Hz band. Boxplots in Figure 3 shows a remarkable separation of the betweenness centrality at right precentral in the low-frequency band, 0.03-0.06Hz, and of the local efficiency at right caudal middle frontal gyrus in the high-frequency band, 0.12-0.25 Hz, between baseline and after receiving the SSRI treatment.
The local efficiency of a brain region is characterized by efficiency of paths from it to other regions, and efficiency of a path is defined by the reciprocal of its length [74]. The nodal clustering coefficient measures the possibility that any two neighbors of the node are also connected. The clustering coefficient is also a measure of functional segregation, which is the ability for specialized processing to occur within densely interconnected groups of brain regions [29]. Reduction in both local efficiency and clustering coefficient at right caudal middle frontal gyrus indicates connectivity inhibition. Participation coefficient and betweenness centrality quantify the hubness of a brain region where the participation coefficient characterizes the diversity of intermodular connections of a node, and betweenness centrality is the fraction (shortest paths passing through the node of interest) of all shortest paths in the network. The raised participation coefficient implies the increase in functional integration of left pars triangularis to the rest of the network. Reduced betweenness centrality at the right precentral gyrus shows dwindling of functional integrity. Within-module degree z-score is a localized, within-module version of degree centrality for regions in a brain network heuristically classified into distinct functional groups. Decreased within-module-degree at the left lateral occipital cortex implies reduced connectivity to other regions in the same functional group.  Based on the finding of local efficiency and clustering coefficient change of right caudal middle frontal gyrus in the 0.12-0.25 Hz network, other regions in the same network are also investigated with a lower criteria (p-value ≤ 0.05). Results, summarized in Table 4, show both decreased clustering coefficient and decreased local efficiency in bilateral caudal middle frontal gyrus, rostral middle frontal gyrus, superior temporal gyrus, left pars triangularis, putamen, and right superior frontal gyrus. No region is found with either increased clustering coefficient or local efficiency. Figure 4 shows the distribution of nodal centrality, degree, of brain regions in all four networks corresponding to the four frequency bands. Node size represents the baseline degree, and node color indicates the direction (increase/decrease) of changes (uncorrected p-value < 0.05) after treatment. The p-values, effect size, and power of the altered regions are listed in Table 5. Results show the decreased degree spread over high-frequency to low-frequency bands in regions: bilateral supramarginal gyrus, bilateral lateral occipital cortex, bilateral superior parietal cortex, left inferior parietal cortex, left superior temporal gyrus, right caudal anterior cingulate cortex, right pars opercularis, right caudal middle frontal gyrus, right precentral gyrus, and right inferior temporal gyrus. Changes in the right hemisphere distribute in anterior cortical regions for the high-frequency (0.03 to 0.25 Hz) bands, and posterior cortical regions in the low-frequency (0.015 to 0.03 Hz) band. Regions found to have increased degree are the bilateral amygdala, bilateral paracentral lobule, left pericalcarine cortex, left middle temporal gyrus, left medial orbital frontal cortex, right lingual gyrus, right accumbens, and right isthmus-cingulate cortex. These regions are all associated with emotion regulation, as reported in the previous studies.  . Treatment-induced changes in nodal degree for brain regions are presented. In the figure, before-treatment nodal degree is represented in node size. Significance and direction of change are presented in node color: red nodes for brain regions with (uncorrected p-value < 0.05) increased nodal degree after treatment, blue nodes for regions with decreased degree, and green nodes for regions without changes. Statistics for changes of red and blue regions are summarized in Table 5 and full name of brain regions are listed in Table 2.

Clinical Correlations
To further understand the functional significance of the significant results listed in Table 3, we conduct correlation analyses with clinical measures. We first investigate if topological measures correlate with BDI scores, which would have suggested relevance to the disease process. The module degree z-score of left lateral occipital cortex in the 0.03-0.06 Hz band show correlation with the BDI scores (both baseline and after treatment) with p-value = 0.02, as shown in Figure 5. Second, we test if the pre-treatment topological measures correlate with the BDI scores and pre-post treatment changes in BDI-II scores to assess if these could serve as treatment outcome predictors. However, no significant correlations were found. Finally, we test correlations between changes in topological measurements and pre-post treatment changes in BDI-II scores, which would have suggested this change as a treatment mechanism. However, no significant correlations were found. Figure 5. Network measure, module degree z-score of left lateral occipital cortex in the 0.03-0.06 Hz band, correlates with clinical BDI score for both baseline and after treatment with coefficient R = 24.4683 and its p-value = 0.020227. Note that the dependence between baseline and posttreatment measurements could cause the inflation of false positive rate. One subject who could not finish the after-treatment clinical assessment is excluded from the regression.

Treatment-Related Changes in Functional Connectivity in Four Frequency Bands
The non-parametric sign permutation test identified (p-value < 0.05 with Bonferroni correction for 2850 connections) decreased connectivity between right fusiform gyrus and left superior temporal gyrus, between right rostral anterior cingulate cortex and right pars opercularis, between right superior temporal gyrus and right medial orbitofrontal cortex, and between right superior parietal cortex and left caudal anterior cingulate cortex in the 0.12-0.25 Hz network; between right medial orbitofrontal cortex and right supramarginal gyrus and between right hippocampus and right rostral middle frontal gyrus in the 0.06-0.12 Hz network; between right thalamus and right superior parietal cortex, between right cuneus cortex and right precentral gyrus, and between right lateral occipital cortex and right pars opercularis in the 0.03-0.06 Hz network; and between left caudal anterior cingulate cortex and right supramarginal gyrus, between left inferior temporal gyrus and left pars opercularis, between left pars opercularis and right fusiform gyrus, and be-tween right caudal middle frontal gyrus and right pars triangularis in the 0.015-0.03 Hz network. The corresponding p-value, effect size, and power are summarized in Table 6, and the boxplots for the distribution of connectivity and its changes are shown in Figure 6. The results show that all the changes involved decreases in connectivity following the treatment. This implies desynchronization of brain activities at the end regions of the listed connections. While considering loosening the correction to accept p-value < 0.0005, eight additional connections are included with decreased connectivity after treatment in Table 6. Notably, most of these connections link regions that are known to be involved in emotional processing and regulation, including hippocampus, thalamus, pars opercularis, fusiform, and various cortical regions in the frontal and temporal lobes, respectively. Moreover, the identified connections are mainly part of the ascending and descending serotonergic pathways, which the SSRIs may affect [75].

Subnetwork Analysis
In addition to identifying changes with Bonferroni correction, a false discovery rate (FDR) approach is used as an alternative to identify altered subnetworks. FDR approach controls for a low proportion of false positives in comparison to the Bonferroni correction, which controls the family wise error rate for the probability of making a false positive conclusion. Figure 7 shows the degree distributions, represented by the size of nodes, for significantly changed subnetworks corresponding to the four frequency bands. The subnetworks are formed by the significantly altered connections determined using the Benjamini-Hochberg procedure for FDR control with α = 0.05. There are 80, 42, 210, and 259 connections in the 0.12-0.25 Hz subnetwork, 0.06-0.12 Hz subnetwork, 0.03-0.06 Hz subnetwork, and 0.015-0.03 Hz subnetwork, respectively. Nodes in the subnetworks include regions with significant topological or degree change presented in Section 3.1, such as the amygdala, caudal middle frontal gyrus, pars triangularis, precentral gyrus, and lateral occipital cortex. Many of the frontal regions including bilateral precentral gyrus, bilateral caudal middle frontal gyrus, bilateral caudal anterior cingulate cortex, bilateral pars opercularis, bilateral rostral middle frontal gyrus, bilateral pars orbitalis, and bilateral lateral orbital frontal cortex are shown associated with numerous significant connectivity changes in the two low-frequency bands. In the two high-frequency bands, changes associated with these regions are more significant in the right hemisphere than in the left hemisphere. Left bank superior temporal sulcus shows significant change across all four frequency bands, and bilateral changes of the superior parietal cortex can be found in the two low-frequency bands. Moreover, the degree of bilateral accumbens increases as the frequency decreases. Additionally, Figure 7b shows minimal changes, Figure 7a shows some minor impacted (medium node size) regions in frontal and temporal areas in the right hemisphere, and they become most significant with the presence of parietal regions when frequency decreases as shown in Figure 7c,d. In contrast, the degree (node size) of regions in the right temporal lobe increases when the frequency decreases.

Clinical Correlation: Baseline Connectivity Predicts Treatment Changes
In Figure 8, the change in clinical assessment (BDI score) is to found correlate with the baseline functional connectivity between right medial orbitalfrontal cortex and right supramarginal gyrus in the 0.06-0.12 Hz band. In other words, higher connectivity in this circuit at baseline was associated with greater clinical improvement. Figure 8. Scatter plot shows correlation between improvement in BDI score and baseline functional connectivity on connection between right medial orbitalfrontal cortex and right supramarginal gyrus in the 0.06-0.12 Hz band (coefficient R = −82.2992 and its p-value = 0.0086528). Age, gender, and IQ were partialed out. One subject who could not finish the after-treatment clinical assessment is excluded from the regression.

Treatment Impact on Network Topology and Functional Connectivity
This study investigates SSRI-induced changes in brain resting-state functional network topology and connectivity using graph-theory-based network analysis. As hypothesized, patients showed changes after treatment in both topological structures and functional connectivity of brain networks. Additionally, the finding that before-treatment connectivity between right medial orbital frontal cortex and right supramarginal gyrus was associated with greater clinical improvement may suggest a promising avenue for research investigating brain predictors of treatment response that could eventually lead to new clinical tools to help clinicians prescribe the most effective medicine for individual patients.
The graph-theory-based network analysis offers a mathematical schematic for characterizing and quantifying the global and regional topological structure of a brain network. To the best of our knowledge, no prior work has investigated treatment-related changes in brain network structure depending on frequencies for MDD. The most relevant work we can find is regarding the disruption of brain network regardless of frequency [46,[76][77][78][79]. Furthermore, to our best knowledge, no prior work has investigated SSRIs-induced topological change in the brain network.

Global Topology
Changes in the global measure of normalized clustering coefficient and smallworldness after treatment are identified. Increased clustering coefficient and smallworldness imply the construction of cliquishness in the network structure. Previous studies have reported decreased clustering coefficient and decreased smallworldness in the MDD patients versus healthy control, which indicates that the treatment change is in the correct direction for improvement [80,81].

Regional Changes and Correlation to Literature
Although thechanges in global measures in Figure 2 are not statistically significant, the changes in local topology support the observed global alterations. These local changes have not yet been reported in any prior SSRI-related study. Therefore, in this section, the changes will be related to previous findings in MDD (may or may not specifically for adolescent).
The reduction in the local efficiency and clustering coefficient at the right caudal middle frontal gyrus in the 0.12-0.25 Hz band is the most significant according to the p-value and effect size in Table 3. The finding is new and, to our best knowledge, no literature has reported MDD-related change in this region for adolescents. The change of the betweenness centrality at the right precentral gyrus in the 0.03-0.06 Hz band is the largest according to the effect size in the table. In the region, in MDD patients, previous studies have reported decreased activity, reduced task-related activation, and increased functional connectivity [82][83][84]. For left pars triangularis, previous studies have reported decreased cortical thickness in depression patients, diminishing pars triangularis in functional networks for depression patients, functional connectivity between pars triangularis and frontal eye field as predictors of acute depression treatment outcome, and pars triangularis as a discriminating biomarker between depression patients and healthy controls [85][86][87][88]. For left lateral occipital cortex, in MDD patients, previous studies have revealed abnormal local intrinsic gray-matter connectivity, decreased baseline blood-oxygen level dependent (BOLD) signal, the correlation between suicidality and connectivity between lateral occipital cortex and fusiform, elevated functional activation, and increased cortical surface [89][90][91][92][93][94]. Even though the previous study has discovered similar correlations in MDD patients and correlation between gray matter thickness in areas of parietal and temporal cortices and antidepressant treatments [95], no finding regarding the treatment-related functional changes for regions in Table 3 has been reported before.
In addition to network efficiency and clustering coefficient, we also investigated the centrality of brain regions using the nodal degree, shown in Table 5, which characterizes the importance of a node in the whole brain network. After treatment, decreased centrality is found at the left inferior parietal cortex, right caudal anterior cingulate cortex, left superior temporal gyrus, bilateral supramarginal gyrus, right pars opercularis, right caudal middle frontal gyrus, right precentral gyrus, bilateral lateral occipital cortex, bilateral superior parietal cortex, and right inferior temporal gyrus. Increased centrality is found at the bilateral amygdala, right lingual gyrus, bilateral paracentral lobule, right accumbens, left pericalcarine cortex, left middle temporal gyrus, right isthmus-cingulate cortex, and left medial orbital frontal cortex. Notably, degree changes are frequency-dependent that some regions tend to have changes in the high-frequency band and the other changes tend to be in the low-frequency band, except for the right amygdala, which has changes across all the frequency bands. Brain regions that show altered clustering coefficient, local efficiency, and nodal degree after treatment include the left amygdala, left medial orbitofrontal cortex, and left paracentral lobule.
In the previous cortical thickness study using the same dataset, the increased thickness of the left medial orbitofrontal cortex and the correlations between clinical improvement and increased cortical thickness of the left lateral orbitofrontal cortex and superior frontal cortex were found [107], which potentially relate to increased connectivity. The bilateral middle temporal gyrus also shows a decreased clustering coefficient and local efficiency but not the nodal centrality. All the regions with the topological change due to the antidepressant are in the frontal-limbic system and the default mode network (DMN), which have been previously implicated in MDD [16,46,[108][109][110][111].
The amygdala is known to play a large part in processing negative emotions, such as the initiation of fear and stress responses [112,113], and is centrally implicated in MDD [21,[114][115][116]. In the amygdala, previous work has reported the volume increase in some studies for patients with depression compared to controls [117], increased BOLD response in unipolar depression [103], and reduced activation after treatment to masked fearful faces [118]. The decreased clustering coefficient, local efficiency, and nodal centrality in the amygdala may reflect an improvement in neurocircuitry. Our previous study using the same dataset but a different approach also discovered changes of resting-state amygdala connectivity, including increased connectivity with the right frontal cortex, right central opercular cortex, and Heschl's gyrus and decreased connectivity with right precuneus, right posterior cingulate cortex, left supplementary motor area and with right precentral gyrus [21,22]. The temporal lobe is implicated in the regulation of emotional states [119,120] and abnormalities in the temporal cortex have been reported in major depressive disorder in several studies [120][121][122][123][124]. For instance, fMRI studies show increased response in the middle temporal gyrus in MDD patients using a paradigm of evoking effect with picture-caption pairs [123] and using negative emotional scripts [119]. The finding in this study shows decreased cliquishness in the clustering coefficient and local efficiency at the middle temporal gyrus in resting-state functional brain networks after treatment.

Frequency Dependency in Regional Changes
From Tables 3 and 4 and Figure 4, locations of altered regions show a trend in frequency: Changes in high-frequency bands (0.12-0.25 Hz and 0.06-0.25 Hz) are mainly located in the frontal and temporal areas, and changes in low-frequency bands (0.015-0.03 Hz and 0.03-0.06 Hz) spread into the parietal and occipital areas. This phenomenon is more obvious in the result of FDR subnetworks showing in Figure 7.
In the FDR subnetworks, regions commonly affected across frequency bands are left banks superior temporal sulcus, right pars opercularis, right precentral gyrus, right pars orbitalis, right lateral orbital frontal cortex, bilateral paracentral lobule, and bilateral superior parietal cortex. Changes at regions including bilateral superior temporal gyrus, fusiform gyrus, inferior temporal gyrus, accumbens, supramarginal gyrus, pars triangularis, rostral middle frontal gyrus, superior frontal gyrus, cuneus cortex, pericalcarine cortex, lateral occipital cortex, and right banks superior temporal sulcus appear to be frequencydependent. In the high-frequency bands, i.e., 0.06-0.12 Hz and 0.12-0.25 Hz, the affected regions are mostly frontal and temporal areas. In the low-frequency bands, i.e., 0.03-0.06 Hz and 0.06-0.12 Hz, parietal and subcortical regions appear affected in addition to the areas identified in the high-frequency networks.
Previous studies have reported subsets of frequency bands as follows: Frequencies between 0.010 and 0.027 Hz may reflect cortical neuronal activity, frequencies between 0.027 and 0.073 Hz may reflect basal ganglia activity, and frequencies between 0.073 and 0.198 Hz and 0.198 and 0.250 Hz have been associated with physiologic noise and white matter signal, respectively [125][126][127][128][129]. Even though there was a debate about whether resting-state fMRI can actually detect white matter activity [130], recent studies tend to support that white matter activity is detectable with resting-state fMRI [131,132].
The frontal lobe is one of the regions that is most consistently identified as associated with MDD [97,133]. Changes in the frontal cortex across all frequency bands may reflect anatomical changes related to white matter. Altered temporal, parietal, and occipital regions in the 0.015-0.03 Hz band are not only primarily related to cortical activities in sensory-motor processing but also participate in a series of higher cognitive functions. The altered network topology in these regions may correspond to potential neural mechanisms mediating the imbalance of these regions' related functional networks.

Functional Connectivity
In order to expand on our past work where we had previously investigated treatmentrelated changes in resting-state functional connectivity between amygdala and all other voxels in the brain [22], we also examined functional connectivity between brain regions in the networks at a global level. With Bonferroni correction, thirteen connections are identified. All the connections in Table 6 (p-value to < 0.0005) show reduced connectivity, and these connections are mostly between limbic, frontal, and temporal regions, as well as the superior parietal cortex.
Moreover, a significant correlation, shown in Figure 8, was found between improvement in clinical assessment and before treatment connectivity of the connection between right medial orbitofrontal cortex and right supramarginal gyrus in the 0.06-0.12 Hz band. Although these correlation analyses are exploratory and the results are not corrected for multiple comparisons, the preliminary findings suggest that the baseline neurological assessment may be used to predict treatment results. MDD is very heterogeneous in that treatment plans have to be individually customized in order to achieve the best result. By better characterizing the neural underpinnings of facts affecting treatment outcomes, this research may pave the way for designing and optimizing treatment for each individual. Taken together, these findings suggest that the network measures derived from graph theory in this study are clinically meaningful, and they may shed light on the neurological underpinnings for better treatment results.
Lastly, we would like to emphasize that this work should be considered as a first and important exploratory study. Treatment-impacted regions, connections, and networks identified in this work have overlap with those implicated in MDD [134]. However, the findings of the SSRI-induced changes and relationships with the clinical assessment reported in this paper require replication with larger data sets. Once validated, these findings could serve as bases for longitudinal studies to answer important questions, such as the following: (1) how network topology and functional connectivity evolve with different lengths of the SSRI intervention; (2) how brain network status correlates with the clinical outcome; (3) how the initial brain network status relates with the intervention outcome; and (4) what and how factors (frequency, strength, etc.) of the intervention contribute to the brain change and clinical improvement.

Methodological Considerations
For the generation of brain functional network and calculation of the network measures, thresholding is utilized to remove insignificant and weak connections that potentially are spurious and could obscure the structure of the significant and strong portion of the network [29]. The threshold is determined per subject in order to ensure the consistent number of connections across subjects for each network density. Any connectivity with a negative correlation coefficient is not considered during the topological analysis, as suggested in the previous study [29].
In this paper, the topological measure-versus-density curves are compared based on the area under curve to test statistical significance, which averages the measure and serves as a scalar summary of the curve values across densities. However, averaging could reduce statistical significance. Another approach is by performing massive comparisons at every single density, which will have higher significance but is also prone to the high type-I error and high sensitivity to noise. To produce a better comprehensive examination of the entire topological structure of the original weighted connectivity graph, in the future study, AUCs for subdivided density ranges or AUCs calculated based on weighted averaging will be investigated to achieve balance between significance and the control of type-I error.
In a connectomic study, type-I error control is necessary for multiple comparisons in finding the changes from 2850 connections in a brain network. A previous study has suggested that, in a connectomic analysis, ruling out unnecessary candidates for comparison using existing knowledge could result in a smaller type-I error and an improved significance [72]. Therefore, in this study, we employed the traditional false-discovery rate (FDR) control procedure to create a significantly impacted subnetwork. Unlike the more restricted family wise error rate control procedure, the FDR procedure ensures the enclosure of all the true positives and offers a guarantee on the ratio of false-positive to the true positive [135]. By excluding the irrelevant part of the network, the preliminary result shows the significantly affected connections between cortical and subcortical regions as part of the ascending and descending pathways that are expected to respond to SSRIs and highlights the subcortical regions associated to these connections. Future investigation will result in a complete study of the topological structure of the subnetwork.
The parcellation of the brain directly determines the size of the brain network and the number of fMRI signals being averaged while computing the seed-based functional connectivity, i.e., Pearson correlation between the average (representative) signal between a pair of regions. In this paper, the commonly used Desikan-Killiany atlas (82 cortical and subcortical regions) is used to define nodes in the functional brain networks. The regions are relatively large compared to the imaging resolution, which means that the advantage of current advanced imaging technology is not fully utilized. Additionally, brain networks derived using different parcellation schemes may show different topological structures. Therefore, it will be interesting to perform similar studies using various atlases and to compare the results.
Additionally, here, functional connectivity is defined based on Pearson correlation coefficients that only characterizes the linear relationship between time series. Some other types of connectivity measures, such as coherence and mutual information, could be utilized to account for a more complicated relationship, such as timing-lag and predictability, between two time series.
Lastly, recent fMRI studies have reported that resting-state functional brain connectivity is highly dynamic [136,137]. In this study, the brain network is formulated as a static network by considering an average representation for the entire period during the resting-state fMRI acquisition, which is consistent with many of the current functional connectomic studies [37,44,46]. However, investigation of the dynamical changes in functional connectivity and network topology is necessary for future work.

Limitations
The small sample size, N = 12, is the major constraint for this study. Although the statistical power for functional connectivity is acceptable, 0.2968-0.9822 (mean 0.6811, SD 0.1838) in Table 6, the power is very limited for global network measures, 0.2834-0.9962 (mean 0.6820, SD 0.1804) in Tables 3 and 4, and for local measures, 0.2110-0.7240 (mean 0.4443, SD 0.1272) in Table 5. Insufficient power would downgrade the ability of discovering a true effect and increase the chance that a significant result does not imply a true effect [138].
In addition to the small sample size, the lack of placebo control that was scanned at similar intervals prevents us from concluding that the changes observed in neural networks are all necessarily due to the treatment, as opposed to spontaneous changes that would have happened anyway. Therefore, the result should be considered preliminary and needs to be validated on large samples in future works.
The absolute values of effect size (Cohen's d) for the five topology measures listed in Table 3 range from 0.7271 to 1.8910 (mean 1.1646, SD 0.4036), indicating large to very large effects [68,139]. Those for connectivity changes that listed in Table 6 range from 0.5818 to 1.6582 (mean 1.0389, SD 0.2697), which indicate medium to very large effects [68,139]. These effect sizes are comparable to the previously reported effect sizes in 26,841 statistical records from 3801 cognitive neuroscience and psychology papers [140]. With the large power for results in Tables 3 and 6, we can eliminate the concern that an effect size could be overestimated by low statistical power [138]. Therefore, even though preliminary findings in this study are promising and could be significant after suitable validation on a large dataset.

Conclusions
In this paper, graph-theory-based complex network analysis has been applied to investigate SSRI-induced alterations of topological organizations and connectivity in resting-state functional brain networks. In the 0.12-0.25 Hz functional brain network, after treatment, MDD patients show decreased local cliquishness and nodal centrality in the middle frontal and superior temporal regions. Furthermore, subnetworks with significantly altered connections present frequency-dependent changes in functional connectivity. Additionally, the baseline connectivity of connection between the right medial orbitofrontal cortex and right supramarginal gyrus in the 0.06-0.12 Hz band is found to be correlated with clinical improvement, which may signal its potential usefulness in predicting treatment outcome, pending confirmation with future studies. The findings of this work may help in gaining new knowledge into the neural underpinnings of the SSRI treatment effect. However, due to the limitation of small sample sizes and lack of placebo control, the results should be viewed as exploratory and need to be validated on large samples in future works. Future efforts will be directed towards studying functional brain networks constructed with different node sets and connectivity measures, exploring the dynamic network structure across time, and testing the results on a larger sample size.