Large-scale dynamic causal modeling of major depressive disorder based on resting-state functional magnetic resonance imaging.

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.

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 (Greicius et al., 2007; 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 well‐established 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 task‐based 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). In particular, by utilizing Bayesian Model Comparison implemented in the PEB framework, one can test different competing hypotheses to reveal the best model.

To overcome the aforementioned limitations of the previous DCM studies and take full advantage of spectral DCM in handling large‐scale EC networks, we conducted an EC study and identified abnormal causal links as potential biomarkers of MDD based on rs‐fMRI data collected from a relatively large sample size (100 [NC] subjects and 100 individuals with first‐episode drug‐naive [FEDN] MDD, to rule out the possible confounding effect of medication). Specifically, we built a relatively large‐scale DCM model containing 27 regions of interests (ROIs) from four probably involved resting‐state networks (DMN, EXE, SAL, and limbic [LIM] networks) using the spectral DCM approach to examine whether MDD alters brain EC and, if so, whether the altered ECs are mainly within‐network and/or between‐network ECs. To achieve this goal, we first applied network‐based statistics (NBS) to detect significant abnormal EC links in MDD. Based on the NBS results, we then formed specific hypotheses, that is, the test within‐ and/or between‐network EC changes in MDD and tested them by using the Bayesian Model Comparison implemented in the PEB framework. Lastly, we performed an automatic search over the reduced PEB models to validate the model comparing results. Results indicated that MDD altered EC both within and between high‐order functional networks. Specifically, at the individual connection level, MDD was mainly characterized by reduced excitatory EC within the DMN and between the DMN and the SAL. At the network level, the inhibitory influence within the DMN was abnormally elevated in the MDD. Overall, our study offers novel insights into the pathophysiological mechanisms of MDD.

METHODS

Participants

A total of 100 FEDN MDD (for convenience, we call them FEDN from now on) patients and 100 NCs were selected from our MDD database (see later for the selection criteria). The demographic and clinical characteristics of participants are shown in Table 1. The FEDN patients were recruited from the psychological counseling outpatient of the First Affiliated Hospital of Guangzhou University of Chinese Medicine, Guangdong, China from September 2015 to June 2018. After an initial screening using a 17‐item Hamilton Rating Scale for Depression (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 self‐reported 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.

Table 1

Demographic and clinical characteristics of participants

Characteristics	FEDN (N = 100)	NC (N = 100)	t/χ 2	p
Age (years)	29.46 ± 9.34a	29.59 ± 10.33	−0.09	.93b
Gender (F/M)	34/66	41/59	1.05	.31c
Education (years)	12.46 ± 3.22a	12.88 ± 2.77	−0.09	.32b
Duration (months)	8.64 ± 10.86a	NA	NA	NA
HDRS‐17	22.15 ± 3.18a	NA	NA	NA

Abbreviations: FEDN, first‐episode drug‐naive major depressive disorder; HDRS‐17, 17‐item Hamilton Depression Rating Scale; NC, normal control.

Mean ± SD.

The p values were obtained by two‐sample t test.

The p value was obtained by a chi‐square test.

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 fluid‐attenuated inversion recovery images. The rs‐fMRI and three‐dimensional T1‐weighted images (3D‐T1WI) in the database have the following parameters: repetition time (TR) / echo time (TE) = 2,000/30, flip angle = 90°, matrix size = 64 × 64, and slice spacing = 1.0 mm for rs‐fMRI, and slice thickness = 1.0 mm, no slice gap, matrix = 256 × 256, field of view = 256 mm for 3D‐T1WI. Because the large‐sample MDD database construction took 4 years, there inevitably existed a few differences in the rs‐fMRI acquisition protocols, such as the field of view, slice thickness, slice number, and total scanning time; to increase sample size, 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 (http://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‐half‐maximum 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.

Table 2

Names and MNI coordinates of 27 regions of interests included in the dynamic causal modeling

Region name		Coordinates (in mm)	Region name		Coordinates (in mm)
DMN (default mode network)			SAL (salience network)
1	Posterior cingulate cortex/Precuneus (PCC_D)	0 −52 7	15	Dorsal anterior cingulate cortex (dACC_S)	0 21 36
2	Medial prefrontal cortex (mPFC_D)	−1 54 27	16	Left anterior PFC (L_aPFC_S)	−35 45 30
3	Left lateral parietal cortex (L_lPar_D)	−46 −66 30	17	Right anterior PFC (R_aPFC_S)	32 45 30
4	Right lateral parietal cortex (R_lPar_D)	49 −63 33	18	Left insula (L_Insula_S)	−41 3 6
5	Left inferior temporal gyrus (L_IT_D)	−61 −24 −9	19	Right insula (R_Insula_S)	41 3 6
6	Right inferior temporal gyrus (R_IT_D)	58 −24 −9	20	Left lateral parietal cortex (L_lPar_S)	−62 −45 30
7	Medial dorsal thalamus (mdThal_D)	0 −12 9	21	Right lateral parietal cortex (R_lPar_S)	62 −45 30
8	Left posterior cerebellum (L_pCERE_D)	−25 −81 −33	LIM (limbic network)
9	Right posterior cerebellum (R_pCERE_D)	25 −81 −33	22	Left subgenual anterior cingulate cortex (L_sgACC_L)	−4 15 −11
EXE (executive control network)			23	Right subgenual anterior cingulate cortex (R_sgACC_L)	4 15 −11
10	Dorsal medial PFC (dmPFC_E)	0 24 46	24	Left amygdala (L_Amyg_L)	−19 −2 −21
11	Left anterior PFC (L_aPFC_E)	−44 45 0	25	Right amygdala (R_Amyg_L)	19 −2 −21
12	Right anterior PFC (R_aPFC_E)	44 45 0	26	Left ventral hippocampus (L_vHPC_L)	−27 −15 −18
13	Left superior parietal lobule (L_sPar_E)	−50 −51 45	27	Right ventral hippocampus (R_vHPC_L)	27 −15 −18
14	Right superior parietal lobule (R_sPar_E)	50 −51 45

Note: The first letter in region name abbreviations (if available) indicates left or right, the last letter indicates network affiliation (D, DMN, E, EXE, S, SAL, L, LIM).

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 are the hidden neuronal states for multiple brain regions, represent any exogenous inputs that are usually absent in the resting‐state fMRI design, is the observed fMRI response, and and 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 ({) and the hemodynamic response function , 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 cross‐spectral density of neuronal fluctuations (Friston et al., 2014): where and represent the Fourier transform of and respectively, and the parameters and 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 is the Fourier transform of , the system's first‐order Volterra kernel (Razi et al., 2015). The expected cross spectra of 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 two‐tailed two‐sample t test was conducted on each connection (link) in the EC network (first testing for NC < FEDN then NC > FEDN, with a p‐value 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 is the observed fMRI data for subject i, is the subject i′ DCM with parameters , and is the observation noise. The function Γ is a second‐level general linear model (GLM) with parameters plus between‐subject variability . Finally, the second‐level parameters is expressed in term of expected value with residual . 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 and Friston 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). It refers to the process of comparing the evidence of the full GLM model with multiple reduced GLMs that have certain combinations of parameters switched off (by fixing the prior expectation at zero; Zeidman et al., 2019b). Such comparisons are very efficient since the evidence and parameters of reduced PEB models can be derived analytically from the full model using Bayesian Model Reduction (Friston et al., 2016). We tested four competing hypotheses developed based on the NBS results: (a) FEDN alters within‐network connections only; (b) FEDN alters between‐network connections only; (c) FEDN alters both within‐ and between‐network connections; or (d) FEDN alters no connections. As summarized in Table 5, each hypothesis contained multiple reduced models that were defined based on the NBS results. Bayesian Model Comparison then assigned a posterior probability (evidence) to each of the candidate models (i.e., the probability of obtaining a particular model given the data) by computing their respective free energy (Zeidman et al., 2019b).

Table 5

Different nested parametric empirical Bayes models with posterior probability for the diagnosis difference

Hypothesis	Model	Connections (ON)	Probability
Within‐network modulation	2	DMN	0
	3	EXE	0
	4	SAL	0
	5	LIM	0
	6	All within‐network	0
Between‐network modulation	7	DMN‐SAL	0
	8	DMN‐SAL; DMN‐LIM; EXE‐SAL	0
	9	DMN‐SAL; DMN‐LIM; DMN‐EXE; EXE‐SAL	0
	10	DMN‐SAL; DMN‐LIM; DMN‐EXE; EXE‐SAL; SAL‐LIM	0
	11	All internetwork	0
Both within‐ and between‐network modulation	1	All connections ON	.27
	12–13	–	–
	14	Models 2 and 9	.04
	15	Models 2 and 10	.03
	16	Models 2 and 11	.03
	17–33	–	–
	34	Models 6 and 9	.34
	35	Models 6 and 10	.26
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.

Lastly, to confirm the Bayesian Model Comparison results, we applied Bayesian Model Reduction to automatically remove redundant EC links from the full PEB model that did not contribute to the model evidence. Specifically, the algorithm implemented a greedy searching over all the permutations of a small set of parameters (i.e., eight parameters) whose removal produces the smallest reduction (i.e., greatest increase) in model fitting (Friston & Penny, 2011). The procedure was repeated until discarding any parameter starts to decrease model evidence or there were no more new parameters to add. After Bayesian Model Reduction, the best 256 reduced models were combined by using Bayesian Model Averaging to account for uncertainty about the underlying model.

To summarize our analysis pipeline, we first estimated EC among the 27 predefined brain regions by using spectral DCM for all 200 subjects (100 NCs and 100 FEDNs). We then applied NBS to identify potentially abnormal EC links in FEDNs and formulated specific hypotheses to test. Next, we took the estimated EC of the full DCM for each subject to the group level and set up a GLM (PEB model) to estimate the effects of group mean and diagnosis difference on each EC link. By utilizing the Bayesian Model Comparison implemented in the PEB framework, we compared multiple reduced models that encoded different hypotheses for finding out the best model that could explain the effect of diagnosis. Finally, we performed an automatic search over the reduced PEB models to confirm the model comparing results.

Data and code availability

Due to sensitive patient information, the fMRI data is protected by the First Affiliated Hospital of Guangzhou University of Chinese Medicine and thus cannot be publicly available. However, upon direct request, we will be able to share all our code for modeling large‐scale DCM for boosting future research and helping to assess the reproducibility of our findings.

RESULTS

There was no significant difference in age (p = .93, two‐sample t test), gender (p = .31, chi‐square test), or education (p = .32, two‐sample t test) between NC and FEDN. The disease duration was 8.64 ± 10.86 months, and the HDRS‐17 score was 22.15 ± 3.18 for the FEDNs (Table 1).

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,b; note that we converted the estimated self‐connection value to units of Hz by multiplying a default value of −0.5 Hz; Zeidman et al., 2019a). The difference in average EC between the NCs and FEDNs (NC–FEDN) is 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).

Figure 1

Network‐based statistics analysis reveals a cluster of effective connectivity (EC) links that show significant group difference between normal control (NC) and first‐episode drug‐naive (FEDN) major depressive disorder. (a) Average effective connectivity from the NC subjects (N = 100). (b) Average effective connectivity from the FEDN subjects (N = 100). (c) Difference in average effective connectivity between the NC and FEDN groups. (d) T‐value from the one‐tailed two‐sample t test for each individual EC link. The links with significant group difference (p < .05, corrected by network‐based statistics) are highlighted in red boxes

The averaged EC strength at the links with significant group differences between NC and FEDN is shown in Table 3. The significant group differences involved 20 ROIs, most of them located in the DMN and SAL (70%, 14/20). The group differences also involved 31 EC links, most of which lied within the DMN (29%) as well as reciprocally between the DMN and SAL (DMN → SAL: 19.4%; SAL → DMN: 22.6%; Figure 2a). The remaining small percentage of EC links was distributed between the SAL and EXE (SAL → EXE: 9.7%), DMN and LIM (DMN → LIM: 6.5%; LIM → DMN: 9.7%), and within the EXE network (3.2%; Figure 2a). In addition, of the 31 EC links, a majority (77.4%, 24/31) maintained excitatory for both NCs and FEDNs, while 19.4% (6/31) switched from excitatory in NC to inhibitory in FEDN and only one (3.2%, 1/31) remained inhibitory for both groups (Figure 2b). For the EC links that remained excitatory or switched from excitatory to inhibitory, the mean absolute EC dropped more than 50% (excitatory/excitatory: from 0.079 to 0.036; excitatory/inhibitory: from 0.041 to −0.015; Figure 2c). By comparison, for the EC link that maintained inhibitory (inhibitory/inhibitory), the EC increased a lot, from −0.002 to −0.05 (Figure 2c). Here, a more negative EC corresponded to stronger inhibition strength. From the group difference result (excitatory/excitatory, the largest category in Figure 2b), the average excitatory EC between NC and FEDN among the four functional networks (DMN, EXE, SAL, and LIM) is separately displayed in Figure 2d, where the EC of all four networks showed a significant decrease in FEDN compared to NC. The LIM network showed the largest reduction (76.8%), followed by the DMN (64.8%). Overall, the excitatory influence was greatly reduced, while the inhibitory effect was substantially increased in FEDN compared to NC.

Table 3

Average effective connectivity of the significant edges between NC and FEDN

Connection	NC	FEDN	Connection	NC	FEDN
DMN → DMN			SAL → DMN
mPFC → R_pCERE	0.12	0.07	L_aPFC → L_pCERE	0.053	0.001
L_IT → mdThal	0.069	0.012	L_Insula → PCC	0.059	0.019
L_IT → L_pCERE	0.106	0.056	L_Insula → R_IT	0.086	0.041
R_IT → R_lPar	0.088	0.034	R_Insula → R_IT	0.079	0.041
mdThal → mPFC	0.093	0.038	R_lPar → R_IT	0.109	0.063
R_pCERE → L_lPar	0.08	0.032	L_aPFC → mPFC	0.046	−0.014
mdThal → L_IT	0.05	−0.01	L_aPFC → R_pCERE	0.037	−0.006
L_pCERE → R_lPar	0.052	−0.005	SAL → EXE
R_pCERE → R_lPar	−0.002	−0.05	dACC → L_aPFC	0.12	0.073
DMN → SAL			dACC → R_aPFC	0.11	0.065
PCC → L_Insula	0.086	0.032	R_Insula → R_aPFC	0.102	0.056
PCC → R_Insula	0.077	0.035	LIM → DMN
R_IT → L_Insula	0.081	0.042	L_sgACC → L_pCERE	0.034	0.008
L_pCERE → L_aPFC	0.062	0.017	R_sgACC → L_pCERE	0.037	0.01
mdThal → R_lPar	0.038	−0.015	R_Amyg → R_lPar	0.035	0.004
L_lPar → L_aPFC	0.021	−0.041	EXE → EXE
DMN → LIM			R_sPar → R_aPFC	0.134	0.083
R_lPar → R_Amyg	0.046	0.010
L_pCERE → R_sgACC	0.038	0.013

Abbreviations: FEDN, first‐episode drug‐naive; NC, normal control.

Figure 2

Distribution of significant effective connectivity (EC) links in the dynamic causal modeling model and change in EC from normal control (NC) to first‐episode drug‐naive (FEDN). (a) Distribution of significant EC links both within and between the four functional networks (default mode, executive control, salience, and limbic networks). (b) Proportion of the significant EC links that remain excitatory (Exc/Exc) or inhibitory (Inh/Inh) under both NC and FEDN conditions and those that change from excitatory in NC to inhibitory in FEDN (Exc/Inh). (c) Difference in average effective connectivity between NC and FEDN for three different types of connections (Exc/Exc, Exc/Inh, and Inh/Inh). (d) Difference in average (excitatory) effective connectivity between NC and FEDN among the four functional networks. The total excitatory (positive) EC (both incoming and outgoing) for each node (in the significant cluster) is summed and averaged within each network. Error bars denote SE. Abbreviations: DMN, default model network; EXE, executive control network; LIM, limbic network; SAL, salience network

The group differences between NC and FEDN was visualized in two weighted and directed graphs using BrainNet Viewer (Supporting Information Figure S2; Xia, Wang, & He, 2013). 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 Table 4. The nodes in the DMN and SAL networks had higher degrees, such as bilateral posterior cerebellum, right lateral parietal cortex, right inferior temporal gyrus, medial dorsal thalamus, left anterior prefrontal cortex, and left insula. At the network level, the SAL network had more outgoing connections than incoming connections, while the EXE network had more incoming connections than outgoing connections showing group differences (Table 4). This indicated that the SAL network influenced other networks more and the EXE network was predominantly under the influence of other networks. Overall, the NBS analysis demonstrated that the default mode and salience networks may be particularly affected by the pathophysiology of MDD.

Table 4

Node degree (number of connections) in the significant cluster identified by network‐based statistics

Node	Total degree	In degree	Out degree	Node	Total degree	In degree	Out degree
DMN				SAL
L_pCERE	7	4	3	L_aPFC	5	2	3
R_lPar	5	4	1	L_Insula	4	2	2
R_IT	5	3	2	R_Insula	3	1	2
R_pCERE	4	2	2	dACC	2	0	2
mdThal	4	1	3	R_lPar	2	1	1
PCC	3	1	2	Average	3.2	1.2	2
mPFC	3	2	1	LIM
L_IT	3	1	2	R_Amyg	2	1	1
L_lPar	2	1	1	R_sgACC	2	1	1
Average	4	2.11	1.89	L_sgACC	1	0	1
EXE				Average	1.67	0.67	1.0
R_aPFC	3	3	0
L_aPFC	1	1	0
R_sPar	1	0	1
Average	1.67	1.33	0.33

To evaluate potential EC changes at the network level, we calculated the average excitatory or inhibitory EC among the four functional networks (DMN, EXE, SAL, and LIM) for both NC and FEDN groups. Specifically, all excitatory or inhibitory EC links (including the 31 links with group differences) among the four networks were separately averaged for each subject and then averaged across subjects for each group. The average excitatory EC for NC and FEDN is 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).

Figure 3

Average effective connectivity for all nodes and connections among the four networks in the dynamic causal modeling model. (a) Average excitatory effective connectivity among the four networks for the normal control (NC) subjects (N = 100). (b) Average excitatory effective connectivity among the four networks for the first‐episode drug‐naive (FEDN) subjects (N = 100). (c) Average inhibitory effective connectivity among the four networks for the NC subjects (N = 100). (d) Average inhibitory effective connectivity among the four networks for the FEDN subjects (N = 100). The inhibitory connection within the DMN shows significant difference between NC and FEDN (p < .05, with FDR correction). Abbreviations: DMN, default model network; EXE, executive control network; LIM, limbic network; SAL, salience network

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 connections best explained the effect of diagnosis using PEB‐based Bayesian Model Comparison (Friston et al., 2016; Zeidman et al., 2019b). The first group of candidate models consisted of connections within the four functional networks: (a) DMN; (b) EXE; (c) SAL; (d) LIM; and also (e) all within‐network connections (Table 5). The second group of candidate models were comprised of particular pairs of functional networks (including both directions): (a) DMN‐SAL; (b) DMN‐SAL, DMN‐LIM, and EXE‐SAL; (c) DMN‐SAL, DMN‐LIM, EXE‐SAL, and DMN‐EXE; (d) DMN‐SAL, DMN‐LIM, EXE‐SAL, DMN‐EXE, and SAL‐LIM; and (e) all between‐network connections (Table 5). That is, the first model corresponded to the largest internetwork connection type (DMN‐SAL) that showed significant group difference based on NBS (Figure 2a). The second model included all internetwork connection types (DMN‐SAL, DMN‐LIM, and EXE‐SAL) that exhibited significant group difference according to NBS (Figure 2a). Given the important role of DMN and SAL identified by NBS, the third model additionally included DMN‐EXE (from the Model 2) and the fourth model additionally included SAL‐LIM (from Model 3), so that all internetwork connections with DMN and/or SAL were included (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 (Table 5). The fourth group of candidate models embodied the null hypothesis that FEDN did not alter any connections (i.e., the effect of 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 difference on each connection (see Section 2). The posterior parameter estimations of the full PEB model are shown in Supporting Information Figure S3a,b, respectively. We then compared the evidence of the full model with the 35 predefined reduced models using Bayesian Model Comparison implemented in the PEB framework (Friston et al., 2016; Zeidman et al., 2019bb). The joint probabilities of all candidate models are shown in Figure 4a, where each element P indicates the posterior probability of the GLM whose parameters for group mean were decided by model i and parameters for group difference set by Model j (Zeidman et al., 2019b). The GLM with the highest probability (34%) corresponded to Model 1 in term of commonalities and Model 34 with respect to group differences (Figure 4a). The model evidence was better visualized when the effect of group means and group differences were unpacked by summing over the individual rows and columns of the joint probability matrix, respectively, to generate the probability of all candidate models for commonalities (Figure 4b) and diagnosis differences (Figure 4c). Clearly, the fully connected model, Model 1, was the only model (100%) that explained the commonalities across subjects (Figure 4b). By comparison, three different models stood out for the group difference, Model 1, Model 34, and Model 35, with a posterior probability of .27, .34, and .26, respectively (Table 5; Figure 4c). Model 34 contained all within‐network connections as well as four internetwork connection assembles (DMN‐SAL, DMN‐LIM, EXE‐SAL, and DMN‐EXE; i.e., Models 6 and 9), while Model 35 additionally included the SAL‐LIM connections (i.e., Models 6 and 10; Table 5). These results suggested that MDD altered a wide range of network connectivities including both within‐network and between‐network connections. The two connection types that were least affected by MDD included SAL‐LIM and EXE‐LIM, since the inclusion of these two types of connections reduced model evidence (i.e., comparing Model 34 with Model 35 and Model 1; Table 5). This was consistent with the NBS results where the significant cluster did not contain SAL‐LIM and EXE‐LIM connections (Figure 2a). In addition, the next three models with moderate evidence included Models 14, 15, and 16 (Figure 4c), which corresponded to DMN and Model 9, DMN and Model 10, and DMN and Model 11, respectively (Table 5). Thus, DMN was particularly targeted by MDD as no other single network resulted in a higher probability when combined with internetwork connections, in good agreement with NBS findings. We then performed a family‐wise analysis by grouping the 36 × 36 = 1,296 models (36 models of the commonalities across subjects, and 36 models of the diagnosis differences) into four families according to the specific type of connections modulated by MDD, as defined earlier (Family 1: within‐network only; Family 2: internetwork only; Family 3: both within‐ and internetworks; Family 4: no modulation). The evidence of the models in each family was then pooled and the families were compared under the assumption that each family was equally likely (Zeidman et al., 2019b). The pooled probability of different combinations of family models is shown in Figure 4d. Clearly, the family models achieved the highest probability (98.9%) when the parameters for both commonalities and group difference were set according to Family 3 (both within‐ and internetwork connections). Hence, the family‐wise analysis confirmed that MDD altered both within‐network and between‐network connectivities.

Figure 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 (Table 5) in terms of commonalities across subjects and differences between subjects due to diagnosis. The best model is Number 1 for the commonalities and 34 for diagnosis difference, with 34% posterior probability. (b) Posterior probability of different models for the commonalities across subjects (summed over the rows of Panel [a]). (c) Posterior probability of different models for the diagnosis difference (summed over the columns of Panel [a]). (d) Pooled probability for different combinations of family models for commonalities and diagnosis differences. Family 1: group level effects on within‐network only (DMN, EXE, SAL, or LIM); Family 2: group level effects on internetwork only; Family 3: group level effects on both within‐networks and internetworks; Family 4: on group level effects. The best family combination (with the highest pooled probability of 98.9%) is Family 3 for both commonalities and diagnosis differences. Abbreviations: DMN, default model network; EXE, executive control network; LIM, limbic network; SAL, salience network

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; Friston et al., 2016; Zeidman et al., 2019b). It removed 25 and 342 connections (out of 729 connections) for commonalities and group differences, respectively, and the posterior parameter estimates of the best 256 models (using Bayesian Model Average) are shown in Supporting Information Figure S3c,d. When a posterior probability (free energy) threshold of .99 was applied, seven connections were further removed, resulting in 380 connections with strong evidence of existed group difference (Figure 5b). For comparison with the NBS‐based results, we plotted the significant EC links identified by NBS in Figure 5a, with the overlap highlighted by red boxes. A majority of the altered EC links (71%, 22/31) found by NBS were included in the reduced PEB model after automatic search. Similar to the NBS‐based results, most of the EC links were found to be reduced in FEDN compared to NC (Figure 5b). However, the PEB identified significantly more altered EC links than NBS did (380 vs. 31 EC links). In addition, all the 27 ROIs were involved in the altered EC links found by PEB, compared to 20 ROIs identified by NBS. This was consistent with the previous Bayesian Model Comparison analysis where the winning models incorporated most of within‐network and between‐network connections (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.

Figure 5

Comparison between the network‐based statistics (NBS) analysis and the parametric empirical Bayes (PEB) model. (a) Significant cluster (in green) identified by NBS. (b) The effective connectivity links that remain in the PEB model (for group difference) after Bayesian Model Reduction. A total of 349 (out of 729) connections are removed. The common connections between NBS and PEB are highlighted with red boxes. Green edges: normal control (NC) > FEDN; blue edges: NC < FEDN

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 task‐based 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.

A triple network model of depression

It has been hypothesized by Menon (2011) that aberrant interactions and engagement of the DMN, EXE, and SAL underlie the psychopathology of a number of neurological and psychiatric disorders including depression, anxiety, and schizophrenia. This triple network 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 self‐referential 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 involving the subgenual anterior cingulate cortex and thalamus and the length of depressive episode correlated positively with the FC in the subgenual anterior cingulate cortex. Zhu et al. (2012) reported dissociation pattern in the DMN in FEDN patients where the FC was increased in the anterior medial cortical regions and decreased in the posterior medial cortical regions. 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 REST‐meta‐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 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). Moreover, while the excitatory EC was significantly reduced within the DMN at the single connection level, the inhibitory influence at the network level was abnormally increased in the DMN in FEDN. Our results suggested that in FEDN, the DMN was overly suppressed with stronger inhibition among different ROIs. Thus, it is possible that both increased FC and inhibitory EC of the DMN contribute to different depressive symptoms. While increased FC contributes to excessive depressive rumination, stronger inhibitory EC may suppress normal self‐referential mental processes to a detrimental level.

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 the acute stress (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). Using high‐model‐order independent component analysis (ICA) of rs‐fMRI data, Manoliu et al. (2014) reported increased FC between SAL and inferior posterior DMN. It should be noted that the above studies estimated unidirectional FC instead of directed EC. Using DCM of rs‐fMRI, we observed reduced EC between the SAL and DMN network, mainly between the insula and posterior cingulate cortex, insula and inferior temporal cortex and between the anterior prefrontal cortex and cerebellum. Our results are consistent with the triple network model proposed by Menon (2011) that weak engagement of the DMN by salient events leads to altered self‐referential mental activity (i.e., excessive rumination) in MDD.

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 (Greicius et al., 2007; 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), and the FC within the EXE (dorsolateral prefrontal cortex) and between the LIM and DMN (amygdala‐posterior cingulate cortex) increased in MDD adolescents (Peters, Burkhouse, Feldhaus, Langenecker, & Jacobs, 2016).

The main difference is that we found reduced (excitatory) EC within the DMN and between the DMN and other networks, while most FC studies reported increased connectivity in these networks, as discussed above. It should be noted that a few FC studies did observe reduced connectivity within the DMN and between particular functional networks and our EC results also support them. For instance, a recent large‐scale multicenter rs‐fMRI study revealed decreased FC within the DMN in recurrent MDD (Yan et al., 2019). Connolly et al. (2017) reported that depressed adolescents showed reduced connectivity between the LIM and executive network (amygdala‐dorsolateral prefrontal cortex) and between the LIM and DMN (amygdala‐ventromedial prefrontal cortex). Deceased FC between the SAL (insula) and frontolimbic networks have also been reported (Guo et al., 2015).

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 between‐network models, plus 945 (15 × 63) possible combinations of within‐network 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 high‐order 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.

Figure S1 The significant cluster identified by NBS is similar when the age, gender and imaging protocol information is regressed out in the one‐tailed two‐sample Edges of the significant cluster without regressing out age, gender and protocol information. (B) Edges of the significant cluster when the age, gender and protocol information is regressed out in the two‐sample t test.

Figure S2. Directed graphs of effective connectivity for both NC and FEDN groups. Only the effective connectivity links with significant group difference (p < .05, corrected by NBS) are shown. The color of the node indicate the affiliation to each of the four functional networks. The size of the node indicates the number of connections (counting both incoming and outgoing connections) and the thickness of the links indicates within‐group averaged strength of EC (taking absolute values). Purple links are excitatory, and black are inhibitory links. Top row: sagittal view; middle row: coronal view; bottom row: axial view. DMN: default model network; EXE: executive control network; SAL: salience network; LIM: limbic network.

Figure S3. Parametric empirical Bayes (PEB) model output. (A) Effective size of group mean for the full model (containing 729 connections); (B) Effective size of group difference for the full model (containing 729 connections); (C) Effective size of group mean for the reduced model (25 connections were removed). (B) Effective size of group difference for the reduced model (342 connections were removed).

Figure S4. Correlation between the T‐value from the one‐tailed two‐sample t test (between NC and FEDN) for each individual effective connectivity (EC) link. (B) PEB model estimate (effective size) for the difference between NC and FEDN. (C) Correlation between t‐value and PET model estimate. Spearman's rank correlation is performed and r is the correlation coefficient.

Click here for additional data file.