{'fraction_non_alphanumeric': 0.037046449244670335, 'fraction_numerical': 0.01660313426514189, 'mean_word_length': 5.048847139197267, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large‐scale dynamic causal modeling of major depressive disorder based on resting‐state functional magnetic resonance imaging

Abstract Major depressive disorder (MDD) is a serious mental illness characterized by dysfunctional connectivity among distributed brain regions. Previous connectome studies based on functional magnetic resonance imaging (fMRI) have focused primarily on undirected functional connectivity and existing directed effective connectivity (EC) studies concerned mostly task‐based fMRI and incorporated only a few brain regions. To overcome these limitations and understand whether MDD is mediated by within‐network or between‐network connectivities, we applied spectral dynamic causal modeling to estimate EC of a large‐scale network with 27 regions of interests from four distributed functional brain networks (default mode, executive control, salience, and limbic networks), based on large sample‐size resting‐state fMRI consisting of 100 healthy subjects and 100 individuals with first‐episode drug‐naive MDD. We applied a newly developed parametric empirical Bayes (PEB) framework to test specific hypotheses. We showed that MDD altered EC both within and between high‐order functional networks. Specifically, MDD is associated with reduced excitatory connectivity mainly within the default mode network (DMN), and between the default mode and salience networks. In addition, the network‐averaged inhibitory EC within the DMN was found to be significantly elevated in the MDD. The coexistence of the reduced excitatory but increased inhibitory causal connections within the DMNs may underlie disrupted self‐recognition and emotional control in MDD. Overall, this study emphasizes that MDD could be associated with altered causal interactions among high‐order brain functional networks.

Overall, this study emphasizes that MDD could be associated with altered causal interactions among high-order brain functional networks.

K E Y W O R D S
brain networks, drug-naive, dynamic causal modeling, effective connectivity, first-episode, major depressive disorder, parametric empirical Bayes, resting-state fMRI

| INTRODUCTION
Major depressive disorder (MDD) is a common but devastating mood disorder that causes severe personal distress and tremendous cost to society. It is the second leading cause of disability worldwide and affects 4.7% of the global population (Ferrari et al., 2013). Despite decades of extensive research, the etiology and pathophysiology of MDD remain not well understood. Converging neuroimaging studies based on blood-oxygen-level-dependent functional magnetic resonance imaging (fMRI) have revealed disrupted functional connectivity (FC) in the resting-state brain networks including the default mode network (DMN), executive control network (EXE), and salience network (SAL; Brakowski et al., 2017;Dutta, McKie, & Deakin, 2014;Menon, 2011;Mulders, van Eijndhoven, Schene, Beckmann, & Tendolkar, 2015;Zheng et al., 2015). For example, MDD was found to be associated with increased FC within the anterior DMN (Zhu et al., 2012) and between the SAL and anterior DMN (Manoliu et al., 2014), and with decreased connectivity between the DMN and EXE (Manoliu et al., 2014). In addition to these "core" triple networks, altered interactions have also been observed in other FC links such as those in the frontolimbic networks Pezawas et al., 2005;Zhong, Pu, & Yao, 2016). It is now generally accepted that MDD can be characterized as a disorder with dysfunctional connections among various brain regions and networks (Drysdale et al., 2017;Mulders et al., 2015;Yan et al., 2019).
However, most existing fMRI studies of MDD focused on undirected synchronizations (i.e., FC) rather than causal influence (directed connectivity) among the neural populations that give rise to the regional fMRI signals (Friston, 2011). Among many other metrics (such as Granger causality analysis, structure equation modeling, psychophysiological interaction), dynamic causal modeling (DCM) is a wellestablished and popular method allowing for modeling such a causal influence (a.k.a. effective connectivity or EC) among different brain regions (Friston, Harrison, & Penny, 2003), which is well suited for MDD study. As the conventional deterministic DCM applies to taskbased fMRI only, two variations of DCM have been developed to estimate EC from resting-state fMRI (rs-fMRI), that is, stochastic DCM (Li et al., 2011) and spectral DCM (Friston, Kahan, Biswal, & Razi, 2014). Stochastic DCM differs from the conventional deterministic DCM in that it models endogenous or random fluctuations in hidden neuronal and physiological states (Li et al., 2011). Because stochastic DCM estimates both EC and hidden neuronal fluctuations, it is computationally intensive and consequently, can only model a limited number of brain regions (Razi et al., 2017). Spectral DCM, on the other hand, estimates the parameters of cross-spectral density of neuronal fluctuations rather than the time-varying fluctuations in neuronal states (Friston et al., 2014). Because spectral DCM does not estimate fluctuations on hidden states, it is much more stable and computationally efficient (Friston et al., 2014). The computational efficiency of spectral DCM is further increased by using FC as prior constraints, making it ideally suited to model large-scale EC networks (Razi et al., 2017).
The majority of DCM studies on MDD have used task-based fMRI yet there are still a handful of reports with rs-fMRI (Hyett, Breakspear, Friston, Guo, & Parker, 2015;Kandilarova, Stoyanov, Kostianev, & Specht, 2018;Li et al., 2017). Hyett et al. (2015) used stochastic DCM to infer spontaneous interactions among five canonical resting-state networks (DMN, EXE, bilateral anterior insula, and left and right frontoparietal attention networks) from 16 participants with melancholia, 16 with nonmelancholic depression, and 16 healthy individuals. It was observed that melancholia was characterized by reduced EC from the anterior insula to right frontoparietal network compared with nonmelancholic depressive subjects and from the anterior insula to EXE in comparison with normal controls (NCs). Li et al. (2017) examined resting-state EC among four regions of the DMN using spectral DCM with 27 MDD patients and 27 healthy controls. The study revealed decreased EC from the left parietal cortex to other DMN regions in the unmedicated patients. More recently, Kandilarova et al. (2018) employed spectral DCM to study EC difference among eight brain regions on a single brain hemisphere between a group of 20 healthy subjects and 20 medicated patients with either MDD or bipolar disorder. They found that depressive patients showed a significantly reduced EC from the anterior insula to the middle frontal gyrus (i.e., dorsolateral prefrontal cortex) and increased EC from the amygdala to the anterior insula.
Although these rs-fMRI studies were able to characterize abnormal EC in resting-state brain networks in MDD, they included a limited number (<10) of brain regions and relatively small (<30 subjects for each group) sample sizes. Such limitations reduced statistical power and hindered a comprehensive understanding of the large-scale EC changes in MDD. In addition, new statistical methods have recently been developed for network analysis and DCM, which offers more power, improved sensitivity, and hypothesis testing capability in MDD study. Parametric empirical Bayes (PEB) is a two-stage empirical Bayesian model dedicated to DCM studies, which considers both the mean and uncertainty (variance) of the EC estimation to infer group differences. PEB was recently used in DCM studies, showing increased sensitivity and robustness (Friston et al., 2016;Zhou et al., 2018) (HDRS-17) with a total score larger than 18 (Hamilton, 1967), two professional psychologists who have more than 10 years experience separately carried out the MDD diagnosis according to the Diagnostic and Statistical Manual (DSM-5, American Psychiatric Association, 2013).
Only the patients who were diagnosed with FEDN by both of the psychologists would be recruited. The inclusion criteria of FEDN are as follows: (a) aged between 18 and 55 years old, (b) right-handed native Chinese speaker, (c) firstly diagnosed with MDD and had no history of any neurological illness or any other forms of psychiatric disorders, and (d) head motion smaller than 2 mm of translation or 2 of rotation in any direction during the rs-fMRI scan (see Section 2.3). Exclusion criteria included (a) a history of significant medical illness, (b) alcohol abuse (a total score ≥ 8 in Alcohol Use Disorders Identification Test [Saunders, Aasland, Babor, de la Fuente, & Grant, 1993]), and (c) contraindications to MRI scan. A rough illness duration of each FEDN patient was selfreported and recorded. NCs were enrolled locally at the same period of time and were physically and mentally healthy based on their medical history and the Mini-International Neuropsychiatric Interview (Sheehan et al., 1998) as well as the total HDRS-17 score of less than 7. This study was conducted in accordance with the Declaration of Helsinki. All participants provided written informed consent and the study was approved by the local ethics committee.

| Image acquisition
MRI data were acquired using a 3.0-T GE Signa HDxt scanner with an 8-channel head-coil within 3 days of diagnosis. The participants were instructed to close their eyes and refrain from thinking anything particular.
Two radiologists made consensus decisions that all participants were free from visible brain abnormalities or any form of lesions based on thick-slice axial T1-and T2-weighted images, as well as a T2-weighted fluidattenuated inversion recovery images. The rs-fMRI and three-dimensional  datasets with two imaging protocols from the same center were used.
However, their influence on the ROI rs-fMRI time series was negligible.

| Image preprocessing
Image preprocessing was performed using SPM12 (www.fil.ion.ucl.ac. uk/spm) and Data Processing Assistant for Resting-State fMRI (DPARSF) version 2.3 (http://rfmri.org/DPARSF) as detailed previously (Yan & Zang, 2010). For each subject, 180 rs-fMRI volumes were remained after removing the first five volumes. The remaining images were corrected for slice acquisition timing and head motion.
Those with excessive head motion were not enrolled (see the inclusion criteria above). The 3D-T1WI was used to guide rs-fMRI registration by using unified segment and Diffeomorphic Anatomical Registration through Exponentiated Lie Algebra (DARTEL) in SPM12.
The rs-fMRI data were smoothed with a 6-mm full-width-at-halfmaximum Gaussian kernel, and further denoised by regressing out several nuisance signals, including the Friston-24 head motion parameters and signals from cerebrospinal fluid and white matter, before linear detrending and temporal band-pass filtering (0.01-0.08 Hz).

| Spectral DCM
We specified a large DCM model with full connectivity consisting of 27 predefined ROIs (see the detailed ROI definition in Table 2) from four distributed high-order cognitive function-related brain functional networks (DMN, EXE, SAL, and LIM) based on their vital role in MDD neuropathology (Drysdale et al., 2017;Dutta et al., 2014;Menon, 2011). The coordinates of the ROIs from the DMN, EXE, and SAL were adopted from Raichle (2011) and those from the LIM network were taken from Drysdale et al. (2017). Please note that we limited the network scale to only 27 ROIs as they are representative nodes in the four functional networks.
Compared to the previous small-scale DCM studies, an EC network with 27 nodes is already relatively large. Our method can be straightforwardly applied to even larger scale EC networks with more computing resources.
After extracting rs-fMRI time series from these ROIs (each of them was a sphere centering at the above coordinates with a radius of 5 mm), we used spectral DCM (Friston et al., 2014; implemented by the spm_dcm_fmri_csd routine in SPM12) to estimate pairwise EC among the 27 ROIs and constructed a directed and weighted graph (representing an EC network) for each subject. The maximal number of iterations was set to be 256 and all DCMs were converged within 256 iterations.
The theoretical foundation and optimization scheme of spectral DCM have been well documented (Friston et al., 2014;Razi et al., 2017;Razi, Kahan, Rees, & Friston, 2015). Here, we briefly recapitulated the basics of spectral DCM for rs-fMRI built on the concept that neuronal fluctuations are driven by intrinsic activity in the absence of external inputs. The generative model of spectral DCM for rs-fMRI is similar to the state-space model in the conventional DCM for task fMRI (Friston et al., 2003) except adding a stochastic term and removing the modulatory component: where x(t) are the hidden neuronal states for multiple brain regions, u (t) represent any exogenous inputs that are usually absent in the resting-state fMRI design, y(t) is the observed fMRI response, and v(t) and e(t) represent endogenous fluctuations and measurement noises, respectively. The matrix A denotes EC among different brain regions, the matrix C represents the influence of exogenous inputs on neuronal states, and θ and φ are the parameters of the neuronal model ({A, C} & θ) and the hemodynamic response function h(x, φ), respectively.
Instead of inverting the stochastic model in the time domain (i.e., stochastic DCM [Li et al., 2011]), spectral DCM converts the stochastic model to a deterministic model by parameterizing the crossspectral density of neuronal fluctuations (Friston et al., 2014): where V(ω) and E(ω) represent the Fourier transform of v(t) and e(t) respectively, and the parameters α v/e and β v/e characterize the amplitude and exponents of the spectral density of neural fluctuations. The state-space model of spectral DCM can be transformed to the following spectral representation: where K(ω) is the Fourier transform of κ(t), the system's first-order Volterra kernel (Razi et al., 2015). The expected cross spectra of y(t) can then be achieved as: where ψ = {θ, φ, α, β} are unknown parameters of the deterministic model that can be estimated using the standard Variational Laplace procedures (Friston, Mattout, Trujillo-Barreto, Ashburner, & Penny, 2007).

| NBS
After EC was estimated individually for all NC and FEDN subjects, the EC parameters of the two groups were compared using NBS for the directed graph to identify potential abnormal EC links in the FEDNs.
We used NBS as an exploratory method to detect abnormal EC links in MDD and to develop hypotheses that can be tested (using a more tailored method, PEB, see next section) due to its popularity in many exploratory brain network studies (Alexander-Bloch, Raznahan, Bullmore, & Giedd, 2013; Wang et al., 2013;Zalesky, Fornito, & Bullmore, 2010;Zhang et al., 2011). Compared to the conventional statistical analysis methods (e.g., two-sample t test), NBS is desirable when performing mass-univariate statistical tests on every connection in a complex graph (Zalesky et al., 2010). The directed version of NBS is based on the same principles as the conventional NBS in that they both rely on finding the connected components. The main difference is that a directed definition of a connected component is defined for the directed version of the NBS, where a pair of nodes must be connected by at least one path, regardless of the direction of the path. The NBS analysis included the following steps (Zalesky et al., 2010). First, a twotailed two-sample t test was conducted on each connection (link) in the EC network (first testing for NC < FEDN then NC > FEDN, with a pvalue of .025 as a threshold for each direction of comparisons). After such linkwise thresholding, a set of suprathreshold connections were identified. Second, topological clusters (if any) among the suprathreshold connections were identified. In each topological cluster, there was at least one path that connected all the involved nodes (all the brain regions were interconnected). Third, a family-wise error rate (FWER)-corrected p-value for each cluster was calculated using a permutation test (conducted for 5,000 times; cluster-level threshold set at p < .05). This is similar to the AlphaSim correction in the voxel-wise statistical analysis adopted by other toolboxes (e.g., SPM and AFNI [https://afni.nimh.nih.gov/]). This statistical analysis method respects the topology of the links, thus well-balancing Type-I and Type-II errors.

| PEB model
PEB is a between-subject hierarchical and empirical Bayesian model that estimates the effects of group mean and group differences on each EC link and it was specifically proposed for DCM studies (Friston et al., 2016). PEB has been recently used to identify altered EC in a small-scale DCM with nine ROIs (Zhou et al., 2018) and its applicability to test specific hypotheses in large EC networks (>25 ROIs) was evaluated in this study for the first time. The PEB model is characterized by the following equations (Friston et al., 2016): where y i is the observed fMRI data for subject i, is the observation noise. The function Γ (2) is a second-level general linear model (GLM) with parameters θ (2) plus between-subject variability ε (2) . Finally, the second-level parameters θ (2) is expressed in term of expected value η with residual ε (3) . In this study, we entered the ECs estimated by spectral DCM approach into the second-level GLM and specified two covariates (overall group mean and between-group difference) in the design matrix. The parameters of the PEB model are estimated using the standard variational Laplace procedures (for more details about the estimation, see Friston et al., 2007 andFriston et al., 2016).

| Bayesian Model Comparison
After estimating the PEB model parameters for the fully connected DCM, we tested specific hypotheses using Bayesian Model Comparison implemented in the PEB framework (Friston et al., 2016;Zeidman et al., 2019b

| NBS-based group difference in EC
The EC of all subjects was individually estimated using spectral DCM approach (see Section 2). From the fully connected model, the EC among the 27 ROIs was compared between the NCs and FEDNs. The average EC matrices for NC and FEDN groups are shown in Figure 1a, b, respectively. We observed that the pattern of the average EC was quite similar between NC and FEDN with relatively higher EC within each functional network compared with internetwork connectivity. In addition, there was strong self-inhibition of each brain region indicated by the dark blue color along the diagonal (Figure 1a shown in Figure 1c where most EC links were weaker in FEDN than NC. Figure 1d shows the t values from the one-tailed two-sample t test for each EC link, with NBS-derived, significant group differences highlighted in red boxes (p < .05 corrected). All these group differences indicated that FEDN had significantly weaker EC compared to NC. To further evaluate the possible influence of gender, age, and different imaging parameters on the group comparison results, we conducted another two-sample t test with gender, age, and imaging protocol information regressed out; the significant cluster identified by the NBS remained almost the same (Supporting Information Figure S1).
The averaged EC strength at the links with significant group differences between NC and FEDN is shown in Table 3 Overall, the excitatory influence was greatly reduced, while the inhibitory effect was substantially increased in FEDN compared to NC.
The group differences between NC and FEDN was visualized in two weighted and directed graphs using BrainNet Viewer (Supporting Information Figure S2; . The node degrees (in and out degrees, as well as the total degree) of the 20 ROIs involved in the group differences were shown in  Figure 3a,b respectively, while the average inhibitory EC for NC and FEDN is shown in Figure 3c,d, respectively. There was no significant difference between NC and FEDN for the internetwork or intranetwork excitatory EC, while the intranetwork inhibitory EC of the DMN was found to be significantly increased (more negative) in FEDN compared to NC (p < .05, FDR corrected).

| PEB-based Bayesian Model Comparison
The above results from NBS indicated that FEDN was associated with EC changes mainly within the DMN and between DMN and SAL, with a smaller percentage of connections within the EXE, between SAL and EXE, and between DMN and LIM (Figure 2a). Based on these observations, we specified four groups of candidate hypotheses (i.e., nested or reduced models; Table 5) to test which mixture of con-  (Table 5).
That is, the first model corresponded to the largest internetwork con-  (Table 5).
The third group of candidate models consisted of all possible combinations of the models in the first group (within-network) and the second group (between-network), with a total of 25 candidate models ( diagnosis was turned off on all within-network and internetwork connections), resulting in a total of 36 candidate models (Table 5).
After specifying the individual nested (reduced) models, we created a second-level GLM to estimate the effects of group mean and group Hence, the family-wise analysis confirmed that MDD altered both within-network and between-network connectivities.
For completeness, in order to confirm whether the NBS approach would give similar results to Bayesian methods implemented in the PEB framework, as well as to ensure that our particular selection of models did not preclude detection of interesting between-group effects, we additionally performed an automatic search over reduced models. This procedure iteratively removed the GLM parameters that did not contribute to model evidence, using Bayesian Model Reduction (see Section 2;  (Table 5 and Figure 4). Lastly, the group differences revealed by the NBS and those by the PEB were highly correlated (r = .39, p < .001, Supporting Information Figure S4). Particularly, the effective size encoding group difference for the EXE-LIM and SAL-LIM connections was relatively weak (Supporting Information Figure S4b), in line with their absence in the best winning model (Model 34; Table 5).
Overall, the PEB method could be more sensitive than the NBS approach and it generally supported the NBS findings.
T A B L E 4 Node degree (number of connections) in the significant cluster identified by network-based statistics

| DISCUSSION
MDD has been increasingly understood as a mental disorder characterized by disrupted interactions among large-scale functional brain networks (Menon, 2011). However, current connectome studies based on rs-fMRI concentrated mostly on undirected connectivity, while existing directed connectivity analysis based on DCM predominantly used taskbased fMRI and often incorporated very few brain regions and small sample size. Such limitations hamper a thorough understanding of the role of within-network versus between-network connectivity in MDD pathophysiology. Our study represented the first attempt to apply spectral DCM to identify impaired causal interactions in large-scale brain network based on rs-fMRI scanned from a large sample size of healthy and FEDN subjects. Notably, we demonstrated convincingly that MDD altered connectivity both within and between high-order functional networks, and detected abnormal causal influences in FEDN among brain regions from a relatively large-scale brain network, which expanded our knowledge on the neurophysiological mechanisms of MDD. It is not a trivial problem when spectral DCM was applied to a large brain network analysis; we have not only demonstrated the feasibility of applying DCM to a large-scale network-based group comparison but also, for the first time, evaluated the complementary utilities of two different statistical analysis approaches to infer group-level difference. Below, we discussed the importance of our findings in the context of existing knowledge of MDD pathophysiology. regulates crucial self-referent, cognitive, and emotional processes that are impaired in depression (Mulders et al., 2015). Numeral fMRI studies have reported abnormal connectivity patterns among the triple networks in MDD (for review, see Brakowski et al., 2017;Dutta et al., 2014;Mulders et al., 2015). Consistently, using DCM of a large-scale network consisting of the triple network plus the LIM network, we identified impaired causal interactions among the four major functional networks. Specifically, at the single connection level, application of NBS revealed a topological cluster that showed decreased EC in FEDN compared to NCs. In particular, most of the abnormal EC links located within the DMN and between the DMN and SAL networks, suggesting that these two functional systems could be especially targeted by FEDN. At the network level, we observed that the intranetwork inhibitory EC within the DMN was abnormally strengthened in FEDN. Our findings agree well with the previous studies that the DMN and SAL networks play a particularly important role in the pathophysiology of MDD (Shao et al., 2018). The topological cluster also contained aberrant nodes in the EXE and LIM networks including the anterior prefrontal cortex, subgenual anterior cingulate cortex, and amygdala, whose functional abnormality in MDD has well been documented (Dutta et al., 2014). It should be noted that the brain regions in the four functional networks often overlap. For example, the subgenual anterior cingulate cortex and hippocampus in the LIM network may be part of the DMN network and the amygdala could be considered as a component of the SAL network (Mulders et al., 2015). Overall, the results of this study support the triple network model as a dysfunctional core network for MDD, whose within-network and between-network connectivities are systematically impaired.

| Altered EC within the DMN
The DMN is highly activated during rest and passive sensory processing and deactivated in cognitively demanding tasks (Greicius, Krasnow, Reiss, & Menon, 2003). Due to its important role in selfreferential processes, the DMN provides the neural substrate for depressive rumination and receives the most attention in clinical MDD research (Hamilton, Farmer, Fogelman, & Gotlib, 2015). A majority of previous studies have reported increased FC in the DMN related to MDD pathophysiology. As the first study to examine the role of the DMN in medicated MDD patients, Greicius et al. (2007) reported Increased FC in the DMN was also observed in MDD adolescents during an emotion identification task (Ho et al., 2015). In contrast, a recently published large-scale multicenter rs-fMRI study (the RESTmeta-MDD project) found that the FC was decreased within the DMN in recurrent MDD, but not in FEDN (Yan et al., 2019). Here, using spectral DCM, we demonstrated that the excitatory EC within T A B L E 5 Different nested parametric empirical Bayes models with posterior probability for the diagnosis difference No modulation 36 All connections OFF 0 Note: Only models with substantial probability are shown for the hypothesis with both within-and between-network modulation.
the DMN was also reduced in FEDN, consistent with a recent rs-fMRI study that reported reduced EC from the left parietal cortex to other DMN regions in MDD patients who were antidepressant drug-free for at least 3 months (Li et al., 2017). In addition, we observed the EC between the cerebellum and other DMN regions substantially decreased in FEDN, in line with previous studies (Guo et al., 2013;Liu et al., 2012)

| Reduced causal interactions between SAL and DMN
The SAL, consisting mainly of the anterior insula and dorsal anterior cingulate cortex, is activated in response to salient stimuli including F I G U R E 4 Comparison of PEB models in a pre-defined model space. (a) Joint probability of all candidate models. The axes list the 36 candidate models (  (Hermans, Henckens, Joëls, & Fernández, 2014;Seeley et al., 2007). It plays an important role in filtering and integrating interceptive, autonomic, and emotional information that is related to affective processes (Seeley et al., 2007), making it a central network in MDD pathophysiology (Menon, 2011). Several studies indicated an increased FC between the anterior DMN and SAL, including the connections between the pregenual ACC and insula (Horn et al., 2010) and between the orbitofrontal cortex and insula (Avery et al., 2013).

| Decreased EC between the SAL and EXE networks
In this study, the EC between the SAL and EXE networks was also found significantly reduced which involved connections from the dorsal anterior cingulate cortex and insula to the anterior prefrontal cortex. Our results accord with two recent DCM studies that showed reduced EC between the insula and executive network in melancholia subjects (Hyett et al., 2015) and decreased influence from the anterior insula to the dorsolateral prefrontal cortex in patients with depression (Kandilarova et al., 2018). Such findings also fit well with the triple network model where weakened interactions from the SAL to EXE network leads to signaling deficits and inappropriate dorsal attentional systems responses (Menon, 2011).

| Relevance to other network measures
As mentioned earlier, a majority of the MDD rs-fMRI studies focused on undirected connectivities (Brakowski et al., 2017;Dutta et al., 2014;Mulders et al., 2015). Our major conclusion that MDD alters EC both within and between multiple functional networks is consistent with existing FC literature. For example, the FC within the DMN was found to be significantly increased in the depressive subjects Ho et al., 2015). Also, it has been shown that the FC between anterior DMN and SAL positively correlated with depression severity in poststroke depression (Balaev, Orlov, Petrushevsky, & Martynova, 2017) Deceased FC between the SAL (insula) and frontolimbic networks have also been reported .
In addition to static FC, dynamic FC (dFC) and graph theory have also been applied to study dysfunctional connectivity in MDD. By comparing unmedicated MDD patients to control participants on dynamic resting-state FC, Kaiser et al. (2016) showed that individuals with MDD were characterized by decreased dynamic (less variable) FC between medial prefrontal cortex (mPFC) and parahippocampal gyrus within the DMN as well as increased dynamic (more variable) FC between mPFC and insula, the latter of which was related to higher levels of recent rumination. Zhi et al. (2018) applied graph theory to examine the impaired topological organization of dFC in MDD and identified five dynamic functional states, three of which showed significant group differences. They observed that MDD patients spent much more time in a weekly-connected State 2, which involved DMN regions and exhibited significantly decreased harmonic centrality involving the parietal lobule, lingual gyrus, and thalamus. The reduced dFC in DMN reported by the two studies (Kaiser et al., 2016;Zhi et al., 2018) is in agreement with our findings that MDD is characterized by lower EC in DMN. In another graph theory study, Zheng et al. (2015) reported that the connectivity degree was increased for the right anterior insula and decreased between the DMN and EXE in MDD, which supported the triple network model (Menon, 2011).

| Complementary nature between NBS and PEB
Our main objective in this study is to examine whether MDD is mediated by within-or between-network connectivities among high-order functional networks. Such hypothesis testing is particularly suited to address using Bayesian Model Comparison implemented in the newly developed PEB framework (Friston et al., 2016;Zeidman et al., 2019b).
However, in the absence of any prior knowledge of EC abnormality in a large network, it is difficult to form specific hypotheses (models) to test due to a large number of possible combinations of within-network and/or between-network connectivities. For example, in this study, there are 15 possible within-network models and 63 possible betweennetwork models, plus 945 (15 × 63) possible combinations of withinnetwork and between-network models. Under such circumstances, NBS is an ideal exploratory method to detect significantly abnormal EC links that could be used to form specific reduced modes to be tested by the PEB approach. Results showed that the hypotheses (reduced models) based on NBS were reasonable. For example, the three best winning models (Models 1, 34, and 35) included all the three significant between-network connectivities identified by NBS (DMN-SAL, EXE-SAL, and DMN-LIM; Table 5). It also indicated that DMN was particularly targeted by MDD compared to other functional networks (see above), consistent with NBS findings.
In addition, our results suggested that the PEB approach could be more sensitive than the NBS, as more group differences were identified by PEB ( Figure 5). Two reasons may underlie the sensitivity of the PEB method. First, the PEB model incorporates the uncertainty (variance) of estimated connection strengths at the single-subject level to make an inference, which increases its sensitivity and makes it robust to outlier subjects with noisy data (Friston et al., 2016). Second, in a large network, many parameters may work together to produce the observed group effects while individual contributions are more difficult to be confidently identified due to the small effect size with large individual variability. Compared with the PEB method, NBS is able to identify a small set of significant EC links that differ between NC and FEDN. Importantly, most of the significant EC links identified by NBS were included in the PEB model, validating the topological cluster of group differences found by NBS. Moreover, the PEB estimated group differences were found to be highly correlated with that from the NBS, indicating the inherent consistency between the two approaches. However, it should be noted that the increased number of group differences found by PEB might also reduce the specificity (more false positives could be revealed). Taken together, we suggest that NBS is suitable to use as an exploratory approach to detect the group differences with relatively larger effect size in a large-scale network that can be used to form specific testing hypotheses, while PEB has the main advantage of testing specific hypotheses and may be more tailored for DCM studies, as previously tested (Zhou et al., 2018).

| CONCLUSIONS
For the first time, a large-scale resting-state DCM revealed systemic causal connectivity changes in subjects with FEDN MDD. We demonstrated that MDD altered both within-network and between-network effective connectivities. In particular, the excitatory influence among the major highorder functional networks was found to be significantly attenuated in FEDN, and the inhibitory influence within the DMN could be abnormally increased. At the regional level, the reciprocal causal connections between the default mode and the salience network may be particularly targeted by MDD, explaining deteriorated self-recognition and emotional control in MDD. This finding has emphasized the role of altered causal interactions among high-order brain functional networks in MDD.

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.