Introduction

The impact of depressive disorders as leading causes of global disease burden was exacerbated during the COVID-19 pandemic (Liao et al., 2021; Renaud-Charest et al., 2021; Santomauro et al., 2021). Patients with major depressive disorder (MDD) display affective and cognitive impairments (Li et al., 2022). As one of the most prevalent of the severe psychiatric diseases, MDD has been studied in depth. Various etiologies have been proposed, such as monoamine neurotransmitter deficiency, altered neuroplasticity and neurogenesis, and structural and functional brain changes (Al Shweiki et al., 2019). However, no single model can explain all aspects of MDD during serial depressive episodes. Therefore, a focus on first-episode drug-naive (FEDN) MDD may better inform the understanding of the pathophysiology of MDD (Lin et al., 2021; Yanmei et al., 2020).

Novel imaging techniques have facilitated neurophysiological research. Functional magnetic resonance imaging (fMRI) may elucidate the etiologic roles of specific brain connectivity networks in the pathogenesis of MDD. For example, MDD patients have exhibited decreased reward network connectivity within the prefrontal-striatal regions (Gong et al., 2017); hyposensitivity in a fronto-parietal network suggestive of altered attention mechanisms (Zweerings et al., 2019); and reduced activation of the left dorsolateral prefrontal cortex (dlPFC), inferring impaired processing of negative emotional stimuli (Trettin et al., 2022). fMRI has been widely used to study both the spatial functional connectivity (FC) and dynamic characteristics of MDD patients. Altered functional dynamics of multiple cortical regions, and aberrant connections between the default mode network (DMN) and the frontoparietal network have been observed in MDD (Demirtaş et al., 2016; Wohlschläger et al., 2018; Zheng et al., 2018; Zhi et al., 2018). In addition, a few EEG and MRI studies have provided insights into dysfunctional spatio-temporal connectivity in MDD (Holmes & Pizzagalli, 2008; Kong et al., 2021; Sheng et al., 2018). However, most studies have yielded varying results regarding the properties of functional networks in MDD. Small datasets, heterogeneous clinical presentations, and evaluations during both initial and recurrent episodes have confounded the replication of results regarding changes of whole brain networks in patients with MDD. Recently, a few studies obtained significant results about FEDN MDD (Wu et al., 2016; Yang et al., 2018; Zhang et al., 2022). The excessive temporal variations of brain FC reflecting abnormal communications between large-scale brain networks over time were present in FEDN MDD (Long et al., 2020). Interestingly, FEDN MDD was not related to abnormalities in the topological architecture of functional brain networks instead of recurrent MDD (Yang et al., 2021). Meanwhile, biomarkers of brain FC have been identified primarily through group-level statistics, which limit their utility in the diagnosis of MDD (Shi et al., 2021).

Compared to traditional statistical analysis, recent studies of large-scale intrinsic functional networks based on machine learning have also demonstrated that disease-specific patterns may diagnose MDD in individual patients (Kong et al., 2021; Patel et al., 2016; Sen et al., 2019; Shi et al., 2021; Wang et al., 2017). Pattern classification techniques applied to large-scale brain networks can distinguish MDD patients from healthy controls (HC) at the individual level.

Consequent to these successful pattern classification studies, we proposed a spatio-temporal machine learning method to identify altered spatio-temporal state patterns of FEDN MDD patients based on resting-state fMRI. Group information guided independent component analysis (GIG-ICA) (Du & Fan, 2013; Du et al., 2017) can identify subtle alterations of disease and can characterize functional networks reliably (Du et al., 2015, 2016; Zhi et al., 2018). In this study, GIG-ICA was used to estimate subject-specific spatio-temporal states of dynamic FCs in a relatively large sample size. A multivariate pattern classification method was used to identify informative intrinsic states and build support vector machine (SVM) classifiers based on informative spatio-temporal states to distinguish FEDN MDD patients from HC. We hypothesized that altered intrinsic networks are changed with corresponding temporal characteristics in FEDN MDD patients. We examined the predictive power of spatio-temporal brain patterns, which might capture more time-varying information and interactions among intrinsic networks.

Methods

Participants

Resting-state fMRI scans were acquired from a big dataset included 543 HCs and 314 FEDN patients from the International Big-Data Center for Depression Research (IBCDR) database supported by the members of REST-meta-MDD Consortium in China (http://yanlab.psych.ac.cn/IBCDR) based on the following selection criteria: 1) subjects in the age range of 18 to 50 years were included; 2) subjects with excessive head motion (mean FD > 0.5 mm or max FD > 2 mm) were excluded. The details of demographic and clinical information are summarized in Table 1. Full details of dataset and its preprocessing pipeline have been described previously (Yan et al., 2019), and some important details about clinical information in each site are summarized in the supplemental material. In addition, data of two sites were assigned to Dataset V for validation, and the data from remaining sites were incorporated into Dataset M for establishing models.

Table 1 Demographic and clinical data
Full size table

Identification of spatio-temporal states

We computed whole-brain dynamic connectivity matrices for each subject between all pairs of regions based on Power's 264 functional ROIs (Power et al., 2011) using a sliding time window method (window length = 60 s) (Allen et al., 2014; Du et al., 2017). All connections were expressed as z scores using Fisher’s z transformation.

GIG-ICA was adopted to window-direction concatenated dynamic connectivity matrices to estimate group-level independent components using an independent healthy dataset (details in the supplemental material), and to extract subject-specific components of the individuals with group-level components as guidance information (Du & Fan, 2013; Du et al., 2017). The number of components was determined to be 50 for exploring fine details of component structures. Each subject-specific independent component was a 264 × 264 connectivity network designated as a spatial state. The time-varying weights course reflecting the importance of the corresponding component was designated as a temporal state. Finally, the dynamic FC of each subject was characterized by 50 spatial–temporal state patterns.

The elimination of site-related effects such as protocol-specific variability is especially important when combining multi-site data for analysis. Consequently, we utilized ComBat-GAM to harmonize the spatial states of subjects to remove site-related effects (Pomponio et al., 2020). ComBat-GAM harmonization was performed using the python package https://github.com/rpomponio/neuroHarmonize.

Supervised feature selection to identify informative states based on spatio-temporal state patterns

A computational wrapper feature selection method (Jing et al., 2019, 2020) that assessed spatio-temporal state utility with respect to a multi-kernel SVM classification algorithm was used to identify the most discriminative combination of states for distinguishing FEDNs from HCs using Dataset M. Similarity measures of spatial states in the SVM model were defined as Riemannian distance on the Grassmann manifold between subjects computed based on their spatial states (Jing et al., 2019, 2020).

Although the informative spatial states were critical for distinguishing patients and HC in previous studies (Jing et al., 2019, 2020; Li et al., 2017), the temporal evolution of MDD might also alter classification. We used the synchronization matrix as a temporal pattern, which was calculated by partial correlation analysis between all pairs of temporal states based on their corresponding temporal weights course, to quantify the temporal characteristic of multiple states. These synchronization matrices were symmetric and positive definite (SPD) matrices, and the collection of SPD matrices forms a cone-shape of a Riemannian manifold as an SPD manifold. Operations to consider whole temporal patterns performed better on the SPD manifold than on Euclidean geometric structures. The distance between two temporal patterns on an SPD manifold can be measured by a Riemannian distance metric.

With a fusion coefficient to combine spatial and temporal pattern, we can build multi-kernel SVM classifiers for spatio-temporal state patterns using nested cross-validation. To avoid bias of the classification, we applied a tenfold cross-validation procedure to optimize the parameters of kernel SVMs.

The details of the approach were shown in the supplemental material.

On the whole, we applied 10-repeated-hold-out cross-validation to identify the discriminative spatio-temporal states with the forward component selection algorithm. In each hold-out cross-validation, we selected 200 patients and 200 controls randomly as the training dataset for balancing classification, and the remaining 255 individuals as the testing dataset. Then, we applied a nested tenfold cross-validation procedure for distinguishing FEDN patients from HC based on each training dataset. Since different states might be selected in different selection procedures due to inter-subject variability, informative spatio-temporal states of FEDN were identified as those that were selected with higher frequency.

Classification between FEDN and HC based on spatio-temporal state patterns

Similarly, the multivariate pattern classification method using nested cross-validation based on informative spatio-temporal states was adopted to FEDN patients and HCs. Similarity measures were defined on the Grassmann and SPD manifolds between subjects computed from their informative spatio-temporal states. In each hold-out cross-validation, the multi-kernel SVM classifier generated a classification score that was a median value of classification scores of its nested tenfold classifiers, with a positive value to indicate FEDN or a negative value to indicate HC.

Nonparametric permutation tests were adopted to estimate statistical significance of the classification results based on the informative spatio-temporal states. Particularly, multi-kernel SVM classification models using disrupted labeled individuals were built to estimate classification performance. This procedure was repeated 200 times. The informative spatio-temporal states were further validated in terms of their classification performance based on individuals from Dataset V.

Comparisons of spatial and temporal characteristics

FC strength differences of the informative spatial states and synchronization coefficients in informative temporal states between FEDN and HC were quantified by using pseudo-two-sample t-tests (permutation test, age, sex and site as covariates).

Correlation analysis between classification scores and clinical responses

A general linear model with outliers removed at 95% confidence intervals was used to investigate correlations among the classification scores, age at onset, and measures of cognitive function (Hamilton Depression Rating Scale [HAMD] and Hamilton Anxiety Rating Scale [HAMA]) in groups, with age, sex, and site as covariates. We excluded subjects with missing clinical information and error classification.

Results

Informative spatio-temporal states in FEDN MDD identified by pattern classification

Informative spatio-temporal states (STSs) included 3 multi-network-coupled states (feature selection frequency > 0.6, Figs. 1, 2A). STS1 featured inter- or intra- FCs located in the default mode network (DMN) and salience network (SAN). STS2 was a hub network centered on the dlPFC and was connected primarily to the DMN, visual, subcortical, and motor-somatosensory networks. STS3 was a relatively complex coupling network, which principally included the visual, DMN, motor-somatosensory and subcortical networks. In STS3, the connections between the DMN or subcortical regions and visual networks were generally negative, while the connections linking DMN and motor-somatosensory networks were positive. Brain connectivity maps weighted by the network degree of these states were mapped into 10 defined common networks based on previous studies (Fig. 1) (Cole et al., 2013; Evan et al., 2016; Power et al., 2011).

Fig. 1
figure 1

Functional connectivity maps of selected informative spatio-temporal states (STS). Gray lines represent the positive connections, and blue lines represent negative connections

Full size image
Fig. 2
figure 2

(A) Selected frequencies of all states. The top 3 states were informative spatio-temporal states (STS). (B) The receiver operating characteristic curves (ROC) and area under ROCs (AUC) of repeated hold-out cross-validations. (C) Distribution histogram of the classification AUCs of permutation tests and the real AUCs of repetitions of cross-validation results

Full size image

The average accuracy of the multi-kernel SVM classifiers built upon these 3 STSs was 73.3% (specificity 74.0%, sensitivity 69.0%) with an area under the receiver operating characteristic curve (AUC) of 0.80 (Figs. 2B). Nonparametric permutation tests showed that the classification results were statistically significant (p < 0.001) as indicated by the histogram of permuted AUCs (Fig. 2C). Furthermore, the classification accuracy of the independent dataset (Dataset V) was 70.2% (specificity 72.9%, sensitivity 66.2%) with an AUC of 0.74.

Statistical comparisons of spatial and temporal characteristics of FEDN MDD and HC

Pseudo-two-sample t-tests were conducted to compare the differences of functional connectivity (permutation test n = 1000, p = 0.01, Fig. 3) and synchronization coefficients (permutation test n = 10,000) of informative spatio-temporal states. For STS1, patients had stronger connectivities between DMN and SAN (precuneus to anterior cingulate cortex) than HC. Relevant complex changes were found in STS2 and STS3. Strengths of the connectivities in STS2 between dlPFC and DMN (precuneus and cingulate cortex) were increased in patients, with other connectivities to dlPFC decreased. In STS3, connectivities to the visual network were generally weaker in patients compared to HC. However, patients exhibited higher inner-connectivities in the motor-somatosensory network and weaker connectivities between the motor-somatosensory and other networks. For temporal patterns, the synchronization coefficient between STS1 and STS3 increased (p = 0.02) in patients compared to HC.

Fig. 3
figure 3

Differences in functional connectivity of informative spatio-temporal states (STS) between FEDN MDD and HC (p < 0.01). Gray lines represent the increased connectivity compared to HC, and blue lines represent decreased connectivity compared to HC

Full size image

Relationship between classification scores and clinical responses

As shown in Fig. 4, patient classification scores were negatively correlated with HAMD (Pearson correlation: r = -0.31, p < 0.00049, sex, age at onset, and site as covariates) and the age at onset (Pearson correlation: r = -0.22, p < 0.011, sex, and site as covariates). No significant association was observed between classification scores and HAMA. Moreover, age at onset was positively related to HAMD (Pearson correlation: r = 0.19, p < 0.026, sex and site as covariates) and HAMA (Pearson correlation: r = 0.34, p < 0.011, sex and site as covariates). Unfortunately, the corrected comparisons after Bonferroni correction were not significant.

Fig. 4
figure 4

Correlation between the SVM scores and clinical characteristics with sex and site as covariates. Scatter plots show that classification scores are negatively correlated with HAMD (A) and age of onsets (B) in the FEDN MDD patients, and age of onsets are positively correlated with HAMD (C) and HAMA (D)

Full size image

Discussion

In the present study, altered dynamic functional patterns estimated by spatio-temporal brain states in FEDN MDD were explored by disease discrimination in a large dataset. FCs and dynamic synchronizations differed significantly between patients and HC. The classification scores of patients correlated with clinical characteristics. Our results suggest that the predictive power of spatio-temporal brain patterns might capture more time-varying information and interactions among intrinsic networks than standard fMRI. We propose that the assessment of spatio-temporal states might represent an insightful clinical and research tool to distinguish between neuropsychiatric patients and HC.

To understand the pathophysiology of MDD, we focused only on FEDN patients. This study adopted an unbiased design for computing group-level intrinsic states in ICA based on 343 HC who were not included in pattern classification. Head motion artifact was minimized by excluding subjects with large degrees of head motion, and by the use of ICA. Instead of evaluating the entire fMRI scan time data assuming stationary FCs (or networks), we applied dynamic FC matrices to capture non-stationary patterns for a better interpretation of the influence of MDD on large-scale brain networks (Miller et al., 2016; Yu et al., 2015). GIG-ICA was used to extract subject-specific brain states based on dynamic FC matrices with stronger independence and better spatial correspondence across subjects (Du & Fan, 2013; Du et al., 2017).

Multivariate pattern classification identified informative spatio-temporal states that differentiated FEDN MDD patients from HC with an AUC of 0.80. ROCs and non-parametric permutation tests demonstrated that the multi-kernel classifiers performed well in distinguishing patients from HC. Validation based on an independent dataset demonstrated that the informative spatio-temporal states were generalizable.

Informative spatio-temporal states included the DMN with SAN (STS1), the dlPFC-hub network (STS2), and a relatively complex coupling network (visual, DMN, motor-somatosensory and subcortical network [STS3]), consistent with previous findings. All three STSs engaged the DMN, which has been associated with various cognitive functions and is preferentially disrupted in MDD (Hamilton et al., 2015; Li et al., 2021; Sambataro et al., 2014). Potential DMN-pattern-based subtypes of MDD have been explored to parse its clinical heterogeneity (Liang et al., 2020). Prominent functional and structural changes have been found particularly in the DMN and SAN; interactions between the SAN and the DMN may be important for cognitive control (Jilka et al., 2014; Shao et al., 2018). We extracted STS1 by a data-driven method, which demonstrated that the DMN-SAN state was critical for dynamic cognitive control in FEDN MDD. The dlPFC region was a component of the cingulo-opercular network and played a key role in the executive control network. Previous studies have indicated that the cingulo-opercular and executive control networks participate in a broad range of cognitive processes in neuropsychiatric disorders (Becker et al., 2021; Coste & Kleinschmidt, 2016; Culbreth et al., 2021; Wu et al., 2016; Zhao et al., 2019). DlPFC dysfunction has also been associated with MDD (Grimm et al., 2008; Meyer et al., 2019; Webler et al., 2022). Moreover, altered connections between dlPFC and DMN may suggest recovery from MDD (Meyer et al., 2019), and circuits connecting the prefrontal cortex and the basal ganglia via the thalamus may represent neuroanatomical substrates of executive processing (Heyder et al., 2004; Menon & D’Esposito, 2022).

Our results described connections between the dlPFC and other brain networks including the DMN, visual, subcortical, and motor-somatosensory networks. These findings illustrate the essential role of the dlPFC in brain dynamic functional patterns. The last informative state (STS3) comprised a more complex coupling network, which primarily included the visual, DMN, motor-somatosensory and subcortical networks. In contrast to the findings of previous studies, the negative connections between the DMN and visual network were decoupled, suggesting that the two networks are anti-correlated in intrinsic functional states. Apart from the DMN, the FCs of the visual, motor-somatosensory and subcortical networks were also consistent with previous findings in MDD (Chen et al., 2019; Kang et al., 2018; Lu et al., 2020; Luo et al., 2021; Yun & Kim, 2021).

The group comparison of FC maps demonstrated that disrupted connectivity was coupled with informative states. For STS1, patients had stronger connectivities between the precuneus and anterior cingulate cortex than HC. The anterior cingulate cortex and the precuneus have been reported to play important roles in the pathophysiology of MDD, although with poor reproducibility of results (Connolly et al., 2013; Lai, 2018; Zheng et al., 2018). For STS2, strengths of the connectivities between dlPFC and DMN were significantly increased in patients, similar to previous results (Ye et al., 2012; Zhang et al., 2022). Notably, there were also a few altered connections between the dlPFC and other regions (superior temporal gyrus in the ventral attention network, inferior occipital gyrus in the visual network, and bilateral superior temporal gyrus in the auditory network); which demonstrated widespread alterations of dlPFC connectivity in patients’ dynamic functional patterns. For STS3, decreased connections exceeded enhanced pathways in MDD. Similar to the findings of earlier studies (Chen et al., 2019; Lu et al., 2020), connectivities associated with the visual network were generally weaker in FEDN patients than in HCs. These findings indicate decreased inter-network FC of the visual network in MDD that suggests a pathogenic role of visual systems (Chen et al., 2019). Our results also suggest decreased FC between the visual and other networks in MDD, and a consequently increased internal role of visual systems. In addition, patients exhibited higher inner-connectivities in the motor-somatosensory network (Kong et al., 2018), as well as weaker connectivities between the motor-somatosensory and other networks.

While some of our results were consistent with those of previous reports, this study had the obvious advantage of using a large multi-site dataset and a data-driven method. Furthermore, we also investigated the group difference of synchronization coefficients between informative states. STS1 and STS3 showed significant differences, and MDD patients had higher synchronization characteristics than HC. The comparison results indicated that the functional specialization of brain networks is disrupted in MDD, and that the inter-network synchronization of DMN-related complex brain states is more vulnerable in neuropsychiatric disorders.

In the present study, classification scores reflected the affinity of fMRI dynamic patterns to FEDN MDD (positive value) or a healthy state (negative value). Because of their association with clinical characteristics, classification scores might serve as useful indicators of altered spatio-temporal FC in MDD. In particular, FEDN MDD patients with higher classification scores had earlier age at onset, suggesting that younger FEDN patients might suffer more severe spatio-temporal functional alterations, and that older patients might experience a mixture of senile neurologic and psychotic disorders in addition to MDD. Notably, age at onset was positively related to HAMD and HAMA scores. Patients with earlier onset MDD tended to have milder symptoms during initial depressive episodes. As a result, classification scores were negatively related to HAMD, which seemed counter-intuitive, but was nonetheless empirically demonstrated.

The present study had several limitations. Clinician-administered assessment scales were limited to the HAMD and HAMA, and clinical information was missing for some individuals from multiple sites. Moreover, no validated questionnaires were applied to assess the symptoms or mental states of either MDD patients or HC. Thus, the assessment scales were not comprehensive enough to explore the intrinsic relationships between informative states and clinical characteristics. The window length and number of components were set empirically with potential impact. Another limitation is that we compared the spatio-temporal states of HCs and FEDN patients at only one time point. A comparison of dynamic functional states of MDD at serial stages of disease progression would be of further interest.

Conclusions

In conclusion, this study applied a multivariate multi-kernel pattern classification method to investigate large-scale spatio-temporal functional brain networks in FEDN MDD. Large-scale dynamic states that included DMN with SAN, a dlPFC-hub network, and a relatively complex coupling network (DMN, visual, motor-somatosensory and subcortical networks) might be altered in MDD. Our results suggest that aberrant spatio-temporal states are informative and predictive in FEDN MDD, and may capture more time-varying information and interactions among intrinsic networks. The assessment of these spatio-temporal states may represent a powerful tool for exploring neuropsychiatric diseases.

Funding Sources

This work was supported by the Beijing Natural Science Foundation (grant numbers 7214299, grant numbers 4214081, grant numbers 7212141, grant numbers 4214080), the National High Technology Research and Development Program of China (grant numbers 2019YFA0706201), the R&D Program of Beijing Municipal Education Commission (grant numbers KM202011232007, grant numbers KM202011232008, grant numbers KM202211232018), and the National Natural Science Foundation of China (grant numbers 82101566). The funding sources had no involvement in the conduction of this research or the preparation of this article.