Aberrant temporal–spatial complexity of intrinsic fluctuations in major depression

Accumulating evidence suggests that the brain is highly dynamic; thus, investigation of brain dynamics especially in brain connectivity would provide crucial information that stationary functional connectivity could miss. This study investigated temporal expressions of spatial modes within the default mode network (DMN), salience network (SN) and cognitive control network (CCN) using a reliable data-driven co-activation pattern (CAP) analysis in two independent data sets. We found enhanced CAP-to-CAP transitions of the SN in patients with MDD. Results suggested enhanced flexibility of this network in the patients. By contrast, we also found reduced spatial consistency and persistence of the DMN in the patients, indicating reduced variability and stability in individuals with MDD. In addition, the patients were characterized by prominent activation of mPFC. Moreover, further correlation analysis revealed that persistence and transitions of RCCN were associated with the severity of depression. Our findings suggest that functional connectivity in the patients may not be simply attenuated or potentiated, but just alternating faster or slower among more complex patterns. The aberrant temporal–spatial complexity of intrinsic fluctuations reflects functional diaschisis of resting-state networks as characteristic of patients with MDD.


Introduction
Major depressive disorder (MDD) is a prevalent psychiatric disorder characterized by the symptoms of depressed mood, anhedonia, cognitive impairments and disturbed sleep or Kaizhong Zheng, Baojuan Li and Hongbing Lu contributed equally.

3
appetite [1][2][3][4]. Researchers have implicated aberrant restingstate functional intrinsic networks in the pathophysiology of MDD [5,6]. Specifically, functional neuroimaging studies have highlighted the involvement of abnormal functional connectivity within the default mode network (DMN), cognitive control network (CCN) and salience network (SN), etc. [7][8][9][10] in MDD. However, findings from previous studies have been mixed. For example, the majority of studies have reported elevated DMN connectivity in the patients [11,12], while some studies also observed decreased DMN connectivity in MDD [13,14]. These inconsistent findings undermine the potential usefulness of resting-state fMRI as an endophenotype or biomarker for major depressive disorder.
One of the possible contributing factors to the contradictory results may be the non-stationary characteristic of brain connectivity [15][16][17][18]. Accumulating evidence has shown that the patterns of brain connectivity could change within time scales of seconds to minutes [16,19,20]. This kind of nonstationary phenomenon has been consistently reported in unconscious anesthetized macaques, healthy controls, as well as patients with Parkinson's disease [21], autism [22], and depression [23]. Thus, traditional (stationary) functional connectivity analysis which simply reports the averaged functional connectivity across the scanning session overlooks dynamics in brain connectivity and hampers the capacity to capture the depression-specific alterations in brain dynamics.
Some attempts have been made to investigate the temporal characteristics of functional connectivity among fixed predefined regions with the sliding-window algorithm which studies the temporal changes of functional connectivity in a truncated sequence of the entire time series [20,24,25]. Altered variability of functional connectivity within the DMN and the frontoparietal network was observed in adults with depression [26]. Moreover, compared to healthy controls, patients with MDD exhibited similar functional connectivity states, but certain states were expressed for markedly different durations [27].
In addition to functional connectivity among fixed predefined regions, the spatial pattern that defines the intrinsic network is also highly dynamic. Recently, Liu et al. proposed the co-activation pattern (CAP) analysis, which provided a richer characterization of functional connectivity that may increase the sensitivity and specificity of dynamic functional network analysis [28,29]. Most importantly, recent studies have shown that CAP analysis could establish the relationship between spatial patterns and temporal changes within intrinsic brain networks, unlike the approach of the slidingwindow [28,[30][31][32]. So far, an increasing number of studies have begun to focus on exploring the dynamic properties of brain connectivity in patients with MDD using the CAP analysis [33,34]. Kaiser et al. investigated the dynamic functional connectivity (DFC) of medial prefrontal cortical (mPFC) regions using the CAP analysis and found decreased DFC between mPFC and the parahippocampal gyrus within the DMN in the patients with MDD [33]. The authors further investigated dynamics of whole-brain states of spatial co-activation and reported abnormal frontoinsular dynamics which were associated with the severity of symptoms in adolescent depression [34]. These studies have clearly shown that CAP analysis was a promising method to uncover the neural mechanisms of major depression.
In the current study, we aimed to investigate the temporal changes in spatial patterns within three large-scale brain networks including the DMN, SN, and CCN in 158 patients with MDD. We proposed a revised CAP analysis algorithm (spCAP). Compared to the original algorithm in Liu et al., the spCAP is able to automatically determine the optimal number of the CAP maps instead of arbitrarily setting it to 8 [28]. In addition, the spCAP employed a spectral clustering algorithm which could accommodate skewed distribution of the sample distance. This algorithm was then used to extract CAPs of the DMN, SN and CCN in patients with MDD. The features including numbers of CAP maps, occurrence rate, within-cluster similarity, persistence and transitions of CAP maps were then compared between the MDD group and healthy controls. More importantly, these findings were reproducible and robust in another independent data set.

Principal data set
A total of 158 patients with major depressive disorder and 102 healthy controls participated in this study. All subjects were recruited from the Zhumadian Psychiatric Hospital. Patients with depression were clinically diagnosed using structural clinical interview for Diagnostic and Statistical Manual of Mental Disorders-IV (DSM-4). The inclusion criteria for MDD groups were as follows: (1) a Hamilton Depression scale (HAMD, 24 items) score of more than 20; (2) no psychotropic medication for at least 2 weeks (6 weeks for fluoxetine) before inclusion. Exclusion criteria for both groups included current or history of systemic medical condition; brain injury; meeting the criteria of DSM-V alcohol/ drug dependence in the past year, or meeting the criteria of DSM-V alcohol/substance abuse in the past 6 months; pregnant or breastfeeding women; having used anticoagulants (heparin, warfarin, etc.), glucocorticoids or medications for thyroid diseases in the past 3 months; abnormal urine toxicology or thyroid screening results and significant current suicidal ideation or suicidal attempt. Table 1 shows the demographical and clinical characteristics of the principal data set.

Replication data set
A total of 99 subjects were recruited including 43 MDD patients and 56 healthy controls. All the subjects were informed about the study's aims and procedures and signed informed consent. The experiment was conducted in accordance with the requirements of the Ethics Committee in Xijing Hospital. Patients with MDD were recruited from the Department of Psychiatry in Xi-jing Hospital. All 43 MDD patients conformed to the diagnostic criteria of DSM-IV for a current episode of MDD as assessed by two experienced psychiatrists. The inclusion criteria were as follows: a Hamilton Depression Scale (HAMD, 24 items) score of more than 18; a of Hamilton Anxiety Scale (HAMA) score of more than 12. Exclusion criteria included: demyelinating disorder, mild cerebral atrophy, previous history of encephalitis, past mental abnormality, vascular malformation, incomplete HAMD test, left-handedness, imaging contraindications due to artificial teeth, deficient cerebellum in fMRI coverage, shadow in T1 image, or problems in spatial normalization, abnormal T1 segmentation and inappropriate head position. Table 1 shows the demographical and clinical characteristics of the replication data set.

Principal data set
Resting-state fMRI data and T1-weighted data were collected in Zhumadian Psychiatric Hospital using a GE Healthcare 3.0 T MR scanner (Signa HDxt scanning, Milwaukee, WI). The subjects were asked to close their eyes, stay awake and lie still in the scanner without performing any explicit task during the fMRI scanning session. The head axis was aligned with the middle axle of MRI machine and two sponges were used to fix the subjects' head to suppress head movement during the experiment. 180 functional images were collected for each subject with the following parameters: TR = 2000 ms, TE = 30 ms, flip angle = 90°, FOV = 220 × 220 mm, slice thickness = 4 mm, spacing = 0.6 mm, matrix size = 64 × 64, number of slices = 33. The high-resolution 3D brain anatomical images using a T1-weighted BRAVO sequence were also collected with the following parameters: TR = 6.8 ms, TE = 2.5 ms, flip angle = 90°, slice thickness = 1 mm, spacing = 0.0 mm, matrix size = 256 × 256, FOV = 256 × 256 mm, number of slices = 192, turnover time (TI) = 1100 ms.
Resting-state fMRI data were preprocessed using the statistical parametric mapping (SPM) software (https:// www. fil. ion. ucl. ac. uk/ spm/ softw are/ spm12/) and the graph theoretical network analysis (GRETNA) toolbox (https:// www. nitrc. org/ proje cts/ gretna/). The first ten volumes of the functional fMRI data were discarded for magnetic saturation. Then, realignment was performed to correct for head motion between fMRI images at different time points by translation and rotation. The high-resolution structural image was then co-registered with functional images and segmented into gray matter, white matter and cerebrospinal fluid (CSF). The deformation parameters from the structural image to the MNI (Montreal Neurological Institute) template were then used to normalize the resting-state fMRI images into a standard space. Additionally, a Gaussian filter with a half maximum width of 6 mm was used to smooth the functional images. Next, each participant's time series were band-pass filtered in the range of 0.01-0.08 Hz and we regressed out the effects of head motion, white matter and cerebrospinal fluid signals. No significant (t = 0.94, p = 0.35) differences in motion/outliers (framewise displacement, FD) were observed between the two groups (FD of MDD = 0.16 ± 0.11, and FD of HCs = 0.14 ± 0.08, separately).
Finally, the signal of each voxel was lowered and then normalized by its standard deviation using the software of MATLAB (https:// www. mathw orks. com/).

Replication data set
Resting-state fMRI data were collected in Xi-jing Hospital using a GE Discovery MR750 3.0 T MRI system. 210 resting-state fMRI images were collected for each subject with the following parameters: repetition time ( Resting-state fMRI data were preprocessed using the statistical parametric mapping (SPM) software (https:// www. fil. ion. ucl. ac. uk/ spm/ softw are/ spm12/) and the graph theoretical network analysis (GRETNA) toolbox (https:// www. nitrc. org/ proje cts/ gretna/). The first ten volumes of the functional fMRI data were discarded for magnetic saturation. The remaining functional fMRI data had been preprocessed for slice timing, correcting for head movement, normalization, smoothing (6 mm), temporally filtering (0.01-0.08 Hz) and head motion, white matter and CSF regression. No significant (t = 1.61, p = 0.11) differences in motion/outliers (framewise displacement, FD) were observed between the two groups (FD of MDD = 0.08 ± 0.05, and FD of HCs = 0.11 ± 0.08, separately). Finally, the signal of each voxel was lowered and then normalized by its standard deviation using the software of MATLAB (https:// www. mathw orks. com/).

Group independent component analysis (GICA)
We used group independent component analysis of fMRI Toolbox (GIFT, http:// mialab. mrn. org/ softw are/ gift/ index. html, v3.0b) to decompose the fMRI images into spatially independent components. The number of independent components was first estimated using the minimum description length (MDL) criterion and we then decomposed the fMRI images into the estimated number of spatially independent components. Finally, the SN, DMN, left and right CCN were identified according to the resting-state network templates created by Smith [35].

Workflow of spontaneous co-activation pattern analysis
Spontaneous activity in resting-state networks-including DMN, SN and bilateral CCN-were investigated using the co-activation pattern (CAP) analysis [28]. Figure 1 shows the workflow of the spontaneous co-activation pattern analysis. In this study, the spontaneous co-activation pattern analysis is mainly divided into five steps: first, we extracted time courses of region of interest (ROI) (Fig. 1A); second, we computed the optimal threshold to select the supra-threshold frames (Fig. 1B); third, the spectral clustering algorithm was applied to cluster the supra-threshold frames (Fig. 1D); fourth, we determined the optimal number of the clusters to obtain CAP maps; fifth, we compared the temporal properties of CAP maps in both groups. Finally, we further extracted the overall dominant CAP maps (dCAP) within resting-state networks and compared the spatial expressions of dCAP in both groups (Fig. 1E).
Extraction of time courses of ROI For each network, we specified a 6 mm-sphere ROI centered on the peak seed locations for the DMN, SN, LCCN and RCCN according to GICA analysis. Time courses from these four ROIs were extracted (Fig. 1A).
Computing the optimal threshold According to a previous study [28], we computed the optimal threshold. Firstly, the spatial correlation (Pearson's correlation coefficient) between the average of the first n% (threshold) frames of each ROI and the corresponding independent components using GICA was computed, resulting in a distribution of the spatial correlation M according to the different threshold. Then we computed the gradient vector G of M: Finally, we computed the optimal threshold n% based on the variation of G.
Spectral clustering algorithm Figure 1C shows the data distribution of patients with MDD in the principal data set using the principal component analysis (PCA). Importantly, we observed that the distance among frames in the principal data set may exhibit skewed distribution which is consistent with previous study [36]. Since the spectral clustering method is independent of the explicit estimation of data distribution [37], we selected the spectral clustering to cluster and pool the supra-threshold frames. The steps of spectral clustering are as follows: Firstly, a cross-correlation matrix J was obtained by calculating the Pearson correlation coefficient between each pair of fMRI frames. A distance matrix S denoted as 1 − J was then used to represent the distances among the suprathreshold frames.
Next, adjacency matrix W was constructed based on S: where σ denotes the width parameter of the Gaussian kernel function and is set to 0.1. We then defined the degree matrix D which is a diagonal matrix: and constructed the Laplacian matrix L = D −1/2 WD 1/2 . The eigenvalues of L were calculated, and the eigenvalues were then sorted in descending order. We obtained the first k (the number of clusters) eigenvalues and calculated corresponding eigenvectors U = [u 1 , u 2 … u k ]. Finally, the K-means algorithm was applied to cluster these frames according to the U matrix.
Determination of the optimal number of the clusters Moreover, we adopted the Calinski-Harabasz (CH) index to identify the optimal number of the clusters [38]. CH is defined as: where n denotes the total number of frames, k denotes the number of clusters, SS W denotes the overall within-cluster variance and SS B denotes the overall between-cluster variance.
Based on a previous study [28], the optimal number of the clusters is not more than 20. Thus, the spectral clustering approach was repeated with k = 2, …, 20, and the optimal number of the clusters was selected on the basis of the highest Calinski-Harabasz values. Finally, the CAP maps were generated through averaging the frames assigned to the same cluster (Fig. 1D).

The temporal properties of CAP maps
To compare between groups on temporal characteristics on each CAP network state, we analyzed CAP analysis with the pooled MDD and HC samples and calculated the measures of . The data were reshaped to a two-dimensional matrix. Then the two-dimensional matrix was descended dimensions using PCA. The first three eigenvectors were selected and displayed. D The extracted fMRI frames were further classified into five groups using the spectral clustering method. Next, the DMN-CAP of each group was generated by averaging all subjects within groups. E The workflow for extracting the overall dominant CAP maps using DMN-CAPs distinct CAP maps' dominance and transitions including overall dwell time in a CAP map, persistence of a CAP maps and total frequency of CAP-to-CAP transitions for each participant. Overall dwell time is defined as the total proportion of frames of different CAP maps to the total number of frames in a subject and persistence is defined as the total frame-to-frame maintenance in a CAP map. Furthermore, we also compared the total frequency of CAP-to-CAP transitions from one CAP map to the other CAP maps. More details were described in a previous study [34]. Finally, we compared the spatial consistency between distinct CAP maps in both groups. Here, the spatial consistency is defined as the association between an individual cluster-assigned time frame and distinct CAP maps.
The spatial expression of the overall dominant CAP maps To select the most reproducible pattern of CAP maps, the overall dominant CAP maps (dCAP) within resting-state networks were extracted (Fig. 1E). Firstly, we computed the temporal fractions (TF) and distinct CAP maps were then sorted in descending order according to TF. The TF is defined as the proportion of the number of fMRI frames assigned to the same cluster to the total number of fMRI frames. Next, the series of temporal frames averages S m was computed: where SM i is the spatial map of CAP i . Then we calculated the Pearson correlation coefficient r s i between SM i and the overall average of all frames. Finally, to remove miscellaneous CAPs, we selected the overall dominant CAP maps (SM i ) through comparing r s i greater than the fixed threshold, where the fixed threshold was chosen as 0.95 according to a previous study [30].

Statistical test
Wilcoxon rank-sum tests (p < 0.05) were applied to determine whether the temporal properties of CAP maps were significantly different between patients with MDD and healthy controls, after the effects of age, gender, and education were regressed out. The relationship between temporal properties of CAP maps and HAMD scores was assessed for the participants with MDD using the Pearson correlation method.

Spatial patterns of the resting-state networks and the time courses of regions of interest (ROIs)
Thirty independent components were obtained by the GICA analyses. The DMN, bilateral CCN and SN were identified according to the spatial maps provided in Smith et al. [35]. The spatial patterns of these networks are shown in Supplementary Fig. 1A

Temporal expression of the CAP maps within three networks
In this study, we examined the measures of distinct CAP maps' dominance and transitions within DMN, SN, LCCN and RCCN by computing overall dwell time in a CAP map, persistence of a CAP map and total frequency of CAP-to-CAP transitions for analysis with the pooled MDD and HC sample. Figures 2 and 3 depict the spatial patterns and the temporal properties of CAP maps within DMN, SN, LCCN and RCCN in both groups. By CAP analysis, we observed that the DMN pattern was decomposed into two CAP maps ( Fig. 2A), the SN pattern was decomposed into three CAP maps (Fig. 2B), the LCCN pattern was decomposed into two CAP maps (Fig. 3A) and the RCCN pattern was decomposed into two CAP maps (Fig. 3B). Increased total frequency of CAP-to-CAP transitions of CAP 3 within SN (p = 0.028, Fig. 2B) in the patients compared with healthy controls were observed in the principal data set (transitions of MDD = 3.38 ± 2.75, and transitions of HCs = 2.78 ± 2.50, separately). Moreover, decreased spatial consistency of CAP2 within DMN (p = 0.018, Fig. 4) was observed in the patients with MDD compared to healthy controls (spatial consistency of MDD = 0.32 ± 0.16, and spatial consistency of HCs = 0.36 ± 0.16, separately). Finally, we calculated the correlations between the temporal properties of CAP maps and HAMD scores for the participants with MDD. However, no significant associations between the temporal properties of CAP maps and HAMD scores were observed in the principal data set. By examining the spatial patterns of the CAP states and the IC map for the DMN (see Fig. 3 of the supplementary materials), we can clearly see that one advantage of the new CAP analysis is that it can provide detailed dynamic properties of spatial modes of the DMN compared with ICA. The DMN component obtained by the traditional ICA was actually highly dynamic and switches between DMN-CAP1 and DMN-CAP2. We also compared the ICA maps of DMN, SN and CCN between the patients and controls (Supplementary Table 3). Compared with the CAPs method, the classical ICA approach was unable to capture the dynamics of brain functional architecture due to its stationary assumption of brain connectivity. In contrast, the CAP analysis could also furnish detailed information on the dwell time, persistence, and numbers of transitions of the CAP states. (a) The DMN pattern was decomposed into two DMN-CAP maps in the principal data set, while that of replication data set was decomposed into three DMN-CAP maps. Significantly reduced persistence of DMN-CAP 3 was detected in the replication data set. (b) The SN pattern was decomposed into three SN-CAP maps in the principal data set, while that of replication data set was decomposed into two SN-CAP maps. Significantly enhanced transitions of SN-CAP 3 was detected in the principal data set. *p < 0.05

Spatial expression of the overall dominant CAP maps
In addition, we further investigated the spatial expression of the overall dominant CAP maps to select the most reproducible pattern of CAP maps. Specifically, the overall dominant CAP, which is the CAP of the maximum temporal fractions, is able to reflect the most repeated spatial pattern dominating the brain repertoire [30]. Figure 6 depicts the spatial expression of the overall dominant CAP maps within distinct resting-state networks. Apparently, their overall dominant CAP maps within resting-state networks are similar (Supplementary Fig. 4).

Reproducible test
In the replication data set, we also used the GICA analyses to obtain 27 independent components. Similarly, the spatial patterns of these networks are shown in Supplementary  Fig. 1B Table 2).
To test the reproducibility of findings from the principal data set, we further repeated the spontaneous co-activation pattern analysis on the replication data set. Similar observations of the number of clusters were observed in the replication data set: the DMN pattern was decomposed into three CAP maps ( Fig. 2A), the SN pattern was decomposed into two CAP maps (Fig. 2B), the LCCN pattern was decomposed intotwo CAP maps (Fig. 3A) and the RCCN pattern was decomposed into two CAP maps (Fig. 3B).
In addition, similar findings were observed in the results of temporal properties. DMN-CAP was detected in reduced pattern (Fig. 4) of spatial consistency within CAP 3 (spatial consistency of MDD = 0.70 ± 1.90, and spatial consistency of HCs = 0.90 ± 2.58, separately) and SN-CAP was observed in enhanced pattern (Fig. 2B) of CAP-to-CAP transition in the patients with MDD compared to healthy controls (transitions of MDD = 4.39 ± 3.05, and transitions of HCs = 3.95 ± 2.62, separately). Moreover, longer persistence of CAP 3 within DMN (p = 0.0294) in healthy controls compared with patients with MDD (persistence of MDD = 1.84 ± 5.74, and PERSISTENCE of HCs = 4.18 ± 6.35, separately) was observed in the replication data set ( Fig. 2A). Similarly, an enhanced pattern of persistence within DMN (persistence of MDD = 4.90 ± 5.74, and persistence of HCs = 5.00 ± 7.09, separately) in healthy controls compared to patients with MDD was detected in the principal data set ( Fig. 2A). In addition, we also analyzed the spatial expressions of the overall dominant CAP maps within different resting-state networks in both groups. Similar spatial expressions of the overall dominant CAP maps of the three networks were detected in the replication data set ( Supplementary Fig. 2).
Moreover, correlation analysis showed that longer persistence of CAP 1 within RCCN was related with decreased HAMD scores (r = − 0.31, p = 0.04, Fig. 5) in the replication data set (persistence of MDD = 6.26 ± 7.95, and HAMD of MDD = 23.4 ± 3.3, separately). In addition, the association between higher frequency of transitions of CAP 2 within RCCN and enhanced HAMD scores was significant (r = 0.32, p = 0.04, Fig. 5) in the replication data set (transitions of MDD = 4.08 ± 2.40, and HAMD of MDD = 23.4 ± 3.3, separately). However, the correlations did not pass Bonferroni correction.

Discussion
In this study, we investigated how temporal expressions of spatial modes were altered in patients with MDD through the CAPs analysis. We found enhanced CAP-to-CAP transitions of SN in subjects with MDD, which suggested enhanced spatial complexity of intrinsic fluctuations. Furthermore, the spatial consistency and persistence of another core brain network (the DMN) was reduced in the patients, suggesting reduced variability and stability in dynamic configuration of DMN. In addition, persistence and frequency of transitions within RCCN was correlated with the severity of depression.
Our findings of enhanced transitions between these CAP modes in the patients suggested increased flexibility of the SN in MDD subjects. These findings are in line with previous studies using dynamic functional connectivity approach that have reported enhanced dynamic functional connectivity between right insula and right precuneus [39] and significantly greater integration of SN in the patients using using resting-state magnetoencephalography (MEG) [40]. Using a data-driven approach, our study further confirmed that flexibility of the SN was disrupted in patients with MDD as reflected by enhanced transitions of recurring coactivation maps. Our results are also in line with previous studies that have consistently reported the involvement of the SN in the neural mechanisms of MDD [5,41]. Abnormal coupling between nodes of the SN-and implicit fractionation of spatial patterns-may accompany a disconnection of functional diaschisis with associated dysfunction of SN in patients, possibly resulting in mood disorders. Furthermore, previous studies have reported that altered the dynamics within insula-CAP maps within salience network using coactivation pattern analysis approach in the participants with high risk of psychosis [42] and autism spectrum disorder [43]. This study further verified aberrant dynamic alternation within SN in the patients with MDD may be common substrates of psychiatric disorder.
Distinct CAP maps correspond to distinct spatial coactivation patterns which are reflective of synchronization of fMRI signals [28]. The patients appeared to mainly stay in a CAP state with prominent activity of the mPFC, while the controls showed longer persistence of another CAP state with prominent activity of the parietal cortex. In addition, reduced stability of a CAP state with prominent activity of the mPFC was observed in patients compared to healthy controls. These findings were in line with previous studies including our own that have reported increased activity or functional connectivity of the mPFC in patients with  [11]. The mPFC is a core region of the anterior subnetwork of the DMN. This region has been shown to be associated with social cognition involving the monitoring of one's own psychological states and mentalizing about the psychological states of others [44,45]. Thus, more prominent mPFC activation and prolonged expression in MDD subjects may account for maladaptive rumination-the process of repetitively and passively thinking about one's negative feelings, possible causes, and consequences in the patients [46,47]. As a matter of fact, recent studies have started to directly investigate association between the temporal patterns of brain activity and depressive symptomatology in patients with MDD. Goodman and co-authors found that increased depressive symptoms were correlated with enhanced frequency and dwell time of the DMN [48]. In addition, another study reported that elevated dynamic functional connectivity between the regions of the DMN and LCCN was associated with more severe symptomatology in depressed patients [33]. In the current study, we observed that persistence and frequency of transitions within RCCN, which was generally involved in cognition control, were correlated with the severity of depression. These findings suggested that investigation of dynamic properties of CAP states may provide useful information for our understanding of the symptomatology in depressed patients.

Limitations
There are some limitations in this work. In the current study, we were interested in three intrinsic brain networks that have consistently been reported to be involved in MDD. Our study demonstrated deficits in brain dynamics of the DMN, SN and CCN. In future studies, we would extend our analysis to other resting-state brain networks such as the attention network, the motor network, and the visual network. Recent studies showed that the functions supported by these networks were also damped in patients with MDD as well [49][50][51]. However, it is not clear whether brain coactivation patterns would also be impaired in the patients. In addition, the subject in the current study were adults between the ages of 18 and 61 years. Although the results are robust in another independent data set, whether our findings are also reproducible in adolescents with MDD or patients with late-onset MDD remains unknown. Our future work will further test the consistency of our findings in different age groups. Moreover, higher HAMD score (p < 0.0001) was observed in the principal data set compared to the replication data set (HAMD of the principal data set = 31.4 ± 7.3, and HAMD of the replication data set = 23.4 ± 3.3, separately), could also contribute to the failure to replication of the group differences. Another limitation of the current study is that we failed to replicate group differences in CAP metrics.
Although similar pattern of changes in dynamic properties could be seen in the replication data set, the between-group differences are not significant. This issue could be due to the relatively small sample size of the replication data set, the clinical heterogeneity of MDD, difference in severity of depression, etc. Finally, we collected the data with traditional rs-fMRI acquisition which could only provide lower temporal resolution of images and limited quality of evaluations of brain dynamics. Advanced multiband multi-echo imaging acquisition will be adopted in our future studies to improve the robustness and repeatability of the CAP method [52].

Conclusions
In conclusion, these results suggest abnormally temporal and spatial changes of functional networks that standard (stationary) functional connectivity could miss; functional connectivity in the patients may not simply be attenuated or potentiated, but just alternate faster or slower among more complex patterns. An abrupt switch of brain patterns may lead to significant considerable variations of standard functional connectivity, which may be a crucial reason for the contradictory results reported in previous studies. The aberrant spatial complexity of intrinsic fluctuations reflects functional diaschisis of resting-state networks as characteristic of patients with MDD.
Research involving human participants and/or animals Human and human samples were used in the current study.

Consent for publication
Last version of the manuscript was approved by the authors before the submission.