Heterogeneous brain dynamic functional connectivity patterns in first‐episode drug‐naive patients with major depressive disorder

Abstract It remains challenging to identify depression accurately due to its biological heterogeneity. As people suffering from depression are associated with functional brain network alterations, we investigated subtypes of patients with first‐episode drug‐naive (FEDN) depression based on brain network characteristics. This study included data from 91 FEDN patients and 91 matched healthy individuals obtained from the International Big‐Data Center for Depression Research. Twenty large‐scale functional connectivity networks were computed using group information guided independent component analysis. A multivariate unsupervised normative modeling method was used to identify subtypes of FEDN and their associated networks, focusing on individual‐level variability among the patients for quantifying deviations of their brain networks from the normative range. Two patient subtypes were identified with distinctive abnormal functional network patterns, consisting of 10 informative connectivity networks, including the default mode network and frontoparietal network. 16% of patients belonged to subtype I with larger extreme deviations from the normal range and shorter illness duration, while 84% belonged to subtype II with weaker extreme deviations and longer illness duration. Moreover, the structural changes in subtype II patients were more complex than the subtype I patients. Compared with healthy controls, both increased and decreased gray matter (GM) abnormalities were identified in widely distributed brain regions in subtype II patients. In contrast, most abnormalities were decreased GM in subtype I. The informative functional network connectivity patterns gleaned from the imaging data can facilitate the accurate identification of FEDN‐MDD subtypes and their associated neurobiological heterogeneity.


| INTRODUCTION
According to the Diagnostic and Statistical Manual of Mental Disorders, 5th Edition (DSM-5), the diagnosis of major depressive disorder (MDD) is made when a patient has any 5 out of 9 symptoms, resulting in dozens of symptom combinations (Goldberg, 2011). MDD is a heterogeneous syndrome varying significantly in illness progression, treatment response, genetics, and neurobiology. Genetic studies have highlighted MDD heterogeneity in differential polygenic risk profiles (Jermy, Glanville, Coleman, Lewis, & Vassos, 2021). The heterogeneity of MDD patients makes it difficult to develop personalized treatments and identify diagnostic biomarkers. To disentangle the heterogeneity in MDD, various attempts have been made to define subtypes based on clinical characteristics Huibers et al., 2015). However, clinical features-based classification is not robust because of the need for a solid theoretical background and neurobiological underpinnings (Drysdale et al., 2017). And also, the clinical symptoms change over time (Schaakxs et al., 2018). The present study aimed to identify subtypes of MDD based on neuroimaging data using machine learning methods.
Recent studies have identified disruptions in functional connectivity (FC) in specific brain networks in MDD (Brakowski et al., 2017;Rolls et al., 2018). Several brain networks were revealed to mediate depressive symptoms, including the default mode network (DMN) and salience brain network (Guo et al., 2014;Manoliu et al., 2013). However, the existing studies yielded inconsistent results in MDD (Runia et al., 2022). For instance, while an earlier study found that MDD was related to increased FC within the DMN (Qin et al., 2015), opposite results were reported by another study (C. Yan et al., 2019). Studies also suggested that the FC was influenced by the treatment. A longitude study found that abnormal FC for MDD patients disappeared after treatment in the posterior DMN, but persisted in the anterior DMN (Li et al., 2013). A similar study reported that 30% of the discriminative FCs were changed towards normal after antidepressant treatment (Qin et al., 2015). This inconsistency may result from the heterogeneity in the patients included in different studies, such as clinical symptoms, illness duration, and medical conditions. Specific patterns of functional brain networks might characterize distinct groups within MDD. Data-driven approaches to explore the MDD subtypes based on neurobiological traits are promising (Loo, Jonge, Romeijn, Kessler, & Schoevers, 2012). Notably, a data-driven clustering method has shown promising performance in grouping chronic, medicated MDD patients into four neurophysiological subtypes based on functional network measures (Drysdale et al., 2017). Two types of medication-naive MDD significantly differed in symptom profiles and brain functional patterns were also identified in a previous study (Wang et al., 2021). A previous study revealed that abnormalities in brain networks in MDD patients could be blurred or hidden by the heterogeneity of the MDD clinical subgroups and brain plasticity may introduce a recovery effect to the abnormal network patterns seen in patients with a relatively short term of illness but unobservable in patients with long term of the illness (Xu et al., 2021). Moreover, a recent study revealed that no significant differences were observed between the chronic, medicated MDD patients and the first-episode drug-naive MDD patients (FEDN-MDD) though both groups were significantly different from healthy controls in their dynamic functional connectivity patterns (Long et al., 2020).
Though many MDD studies assumed that the brain FC remains stationary during the scanning period, brain FC fluctuates over time could provide potentially new information on the brain function (Demirtaş et al., 2016;Hutchison et al., 2013;Mokhtari, Akhlaghi, Simpson, Wu, & Laurienti, 2019;Tu et al., 2019) and recent studies used sliding time window method found that temporal variability and characteristic temporal path length were significantly correlated with depression severity in MDD patients (Kaiser et al., 2016). Several studies have also found increased temporal fluctuations in the DMN (Demirtaş et al., 2016;Kaiser et al., 2016;Long et al., 2020), and one study reported the opposite result (Zhu et al., 2020). Most of these studies selected regions of interest based on previous studies (Cullen et al., 2009;Lui et al., 2011); in contrast, a whole brain analysis would be more objective and least biased to experience.
On the other hand, it has been shown that conventional group difference methods are not equipped to effectively characterize individual differences within groups . Normative modeling is an alternative yet powerful means for characterizing normal variation of healthy controls and elucidating individual patients' deviation from the controls (Bethlehem et al., 2020;Lv et al., 2020;Marquand, Rezek, Buitelaar, & Beckmann, 2016;Shan et al., 2022;Zabihi et al., 2019). Using normative modeling, researchers found that autism patients displayed highly individualized deviations in their cortical thickness, and the deviations were correlated with the severity of repetitive behaviors . Another brain structure study also revealed high heterogeneity in patients with schizophrenia (Lv et al., 2020).
To elucidate MDD heterogeneity and identify robust MDD subtypes based on their functional magnetic resonance imaging (fMRI) data, we adopted the normative modeling method in conjunction with a data-driven unsupervised method to characterize personalized functional connectivity patterns of FEDN-MDD patients and identify FEDN-MDD subtypes, as illustrated in Figure 1. Specifically, we used a data-driven approach, group information-guided independent component analysis (GIG-ICA) (Du et al., 2017), to compute dynamic func- In this study, we employed data from International Big-Data Center for Depression Research (IBCDR), which was collected in China (The data is available by request from any qualified investigator. Website: http://yanlab.psych.ac.cn/IBCDR) (Deng et al., 2022;Ding et al., 2022). The data use agreement was attached as supplementary F I G U R E 1 Methodological overview. (A) Resting-state fMRI data were collected to create dynamic functional connectivity patterns based on the AAL atlas using a sliding time window method. (B) Healthy individuals from an independent dataset were used to compute the subjectspecific networks at the group level, and the back reconstruction algorithm of GIG-ICA was used to compute inherent independent components (ICs) with the time-varying weights of the remaining participants. (C) For the characteristics of healthy-base controls (blue dots), the normative models were estimated between any characteristics and age to identify the normative distribution. Then, the deviations of each FEDN-MDD individual (red dots) were quantified based on these normative models. (D) The deviations of each FEDN-MDD individual were used as features to implement clustering analysis and identify the informative ICs. (E) Comparisons were performed on the clinical, functional connectivity and gray matters across the groups. material. Deidentified and anonymized data were contributed from studies approved by local Institutional Review Boards. All study participants provided written informed consent at their local institution.
Full details on the datasets and preprocessing pipeline have been described previously (C. Yan et al., 2019). Resting-state fMRI was acquired at each site (TR = 2000 ms, TE = 30 ms, Flip Angle = 90 , number of volumes ≥210, opened eye) and was preprocessed using the DPARSF software with a standardized protocol (http://www. rfmri.org/DPARSF). According to the information provided by IBCDR, we included subjects for statistical analyses through the following criteria: (1) patients with first-episode drug-naive (FEDN) depression disorders; (2) subjects with age in the range of 18 to 50; (3) subjects with less head motion (mean FD < = 0.5 mm or max FD < = 2 mm). To avoid the center effect, 164 healthy individuals and 91 FEDN-MDD patients (60 females) from the same center were included.
To prevent bias from these unmatched data, we used the matchIt function in R (https://www.rdocumentation.org/packages/MatchIt/ versions/4.4.0/topics/matchit) to select the comparable group from the healthy controls as the healthy-base group for the patients on the following characteristic: sex, age, and educational levels. Another independent sample (125 control individuals, 71 females) from the other two centers with similar data acquisition was a healthy-ICA group to create the group-level components as guidance information in ICA. Acquisition parameters of the dataset were shown in the supplementary materials (Table S1).

| Dynamic functional connectivity pattern and feature extraction
For each subject, we computed whole-brain dynamic functional connectivity matrices based on m (m = 90) ROIs from the automated anatomical labeling (AAL) template using a sliding time window method (window length = 60s, step size = 1) (Allen et al., 2014).
GIG-ICA was adopted to the dynamic functional connectivity matrices to extract the subject-specific independent components (ICs) and the corresponding time-varying weights reflecting the variability of the ICs (Du et al., 2017). First, we chose the healthy-ICA group to compute the ICs at the group level. Then these group-level ICs were used as guidance information to calculate subject-specific ICs of the FEDN-MDD and healthy-base controls. The number of IC was empirically determined to be 20. Finally, the dynamic functional connectivity pattern of each subject was represented by 20 IC measures (global efficiency coefficient of each IC) and 20 fluctuation coefficients (average change of the frame-wise variation of each IC's time-varying weights) ( Figure 1B). The details of the approach are shown in the supplemental material.

| Normative modeling and estimation of the feature deviations
After obtaining the 20 IC measures and 20 fluctuation coefficients of each healthy-base control, we built the second-order polynomial regression normative model (see supplementary materials) on the relationship between age and these 40 features, respectively. We estimated the 50 th (median) percentile with bootstrapping confidence intervals to quantify the extent of deviations for any FEDN-MDD patient from the normative model. To avoid bias in the model, we used the range between the 5 th and 95 th percentiles for a given age as a scaling coefficient to normalize the z-scores of the feature deviations. Thus, each FEDN-MDD patient of the same age was positioned and a z-score for quantifying each feature was created from the normative range as follows: For individuals with index i, z i jf was the normalized deviation value for feature f of IC j . C:real was the real value of the feature, and C:5 th , C:50 th and C:95 th were the predictions in 5 th , 50 th and 95 th percentiles for the same age respectively.
To assess generalization, we applied 10-fold cross-validation to

| A multivariate unsupervised method to identify the informative features
To investigate the importance of the features and remove the redundancy, a forward selection technique based on the k-means algorithm was used to select features. The selected features were informative in exploring the homogeneous subgroups in the FEDN-MDD cohort.
Specifically, the forward feature selection procedure was carried out 100 runs repeatedly to avoid any bias in the clustering. Since different features might be selected in different runs due to the random selection of initial centers in the k-means algorithm, the selected informative features of FEDN-MDD were identified as those with higher frequency. The details of the approach, including the optimal number of clusters were shown in the supplemental material. The number of clusters in k-means was eventually identified as k = 2, considering the performance and the limitation of the sample size.
Similarly, a 100-repeated clustering procedure based on the selected informative features was applied to cluster the FEDN-MDD patients, and the overall performance was quantified based on these 100 runs. Finally, each FEDN-MDD patient was given 100 individual clustering labels that yielded 100 repeated runs. A certainty measure to evaluate the clustering reliability for each subject based on 100 repeated runs. Notably, the clustering certainty was measured by , where n i is the number of labels i and n is the total number of repeated runs). A higher value indicated higher clustering reliability and vice versa. We excluded 2 subjects with lower clustering certainty (<0.6) in post hoc analyses.
The FEDN-MDD patients' labels were obtained by i ¼ arg max i¼1, …,k ni n À Á as a robust measure to characterize the heterogeneity of FEDN-MDD, and two subtypes of FEDN-MDD were identified.
Pseudo two sample t-tests were conducted to explore the functional connectivity differences of the ICs between healthy-based controls and subtypes of FEDN-MDD using statistical nonparametric mapping software (SnPM13, http://warwick.ac.uk/snpm; permutation tests n = 5000, age, sex and education as covariates).

| Voxel-based morphometry comparisons
To explore if functional changes were associated with structural alterations, voxel-based morphometry (VBM) analysis was carried out to compare gray matter (GM) maps of FEDN patients and healthy subjects. Pseudo two sample t-tests were adopted to compare GM maps between healthy-base controls and subtypes of FEDN-MDD using permutation tests (n = 10000, age, sex, education and total incranial volume (TIV) as covariates).

| Demographic characteristic
Sample characteristics are provided in Table 1. We used the matching method to match the healthy controls from the patients on the following characteristic: sex, age, and educational levels. The included healthy-based group did not differ in sex, age, and educational levels from the FEDN-MDD.

| Informative features for FEDN-MDD
Twelve features from 10 ICs were selected as informative features to cluster the FEDN-MDD. These selected features primarily consisted of IC measures of IC 1 to IC 9, and fluctuation coefficients of IC 2, IC 6, and IC 10 (see Figure 2B). These informative ICs involved the brain network in the default mode network (DMN) (IC 1), the frontoparietal networks (IC 2), the ACC-hub network (IC 3), the visual network (IC 4), the language-related network (IC 5), the bilateral executive control network (IC 6 and IC 7), the hippocampus-hub network (IC 8), the thalamus-hub network (IC 9), and the caudate-hub network (IC 10, for details, see Figure 2A). As shown in Figure 2C, across the 100 runs, the labels of each patient were consistent, with only two patients with lower label consistency excluded (clustering certainty <0.8). The high consistency reflected the clustering reliability for each subject. As shown in the two-dimensional plot in Figure 2D, the plot contained 2 clusters arranged in the first two dimensions.

| Clustering analysis and clinical symptom comparisons in two subtypes of FEDN-MDD
The proposed multivariate unsupervised approach identified two subtypes of the FEDN-MDD, with 14 patients as subtype I and 75 patients as subtype II. The mean silhouette index (Rousseeuw, 1987) of the clustering was 0.461, and the instability of the cluster centers (De Mulder, 2014) was 4.08. The mean stability of the clustering was 0.97. There was a clear distinction between the identified two subtypes on the 12 informative features (see Figure 3A). The average deviations of each feature for each subtype were shown in Figure 3B, and two subtypes showed significant differences in deviation characteristics, especially the 20 IC measures' deviations. To some extent, the deviations pattern of subtype I was more The comparison was conducted on these two groups.
significant than subtype II, and this pattern was limited to the informative IC measures' deviations. A comparison of intra-IC functional connectivity of 10 informative ICs showed a significant difference in each brain network among the three groups (5000 permutation tests, p < .01). The differences are shown in Figure 4. There were more abnormalities across all informative networks in the subtype I patients The severity of depression and anxiety did not differ between subtypes. However, the illness duration in subtype II (22.1 ± 23.92) patients was longer than in subtype I (11.75 ± 12.09) (p = .03, Table 2). In addition, the subtype II patients had fewer extreme deviations in the IC measures and fluctuation coefficients from the normal range.

| Voxel-based morphometry comparisons
Considering the well-documented relationship between the illness duration and the brain structural changes (Kanaan et al., 2009), we also analyzed the differences in whole-brain GM between the two subtypes. Compared to healthy controls, the subtype II patients showed increased GM (mainly in the middle frontal regions and several posterior brain regions) and decreased GM (mainly in temporal and superior frontal regions) distributed widely in the brain. In contrast, subtype I patients, mainly showed decreased GM in the temporal cortex and anterior cingulate gyrus (see Figure 5).

| DISCUSSION
The present study suggested that FEDN-MDD patients might be a heterogeneous group, mainly classified into two subtypes by a multivariate unsupervised method based on normative models. Two subtypes of patients showed different illness durations, GM abnormalities, and functional connectivity alterations in several informative connectivity networks. The subtype I patients were those with more extreme deviations in the functional characteristics of brain networks, and also with shorter illness durations. Subtype II patients had F I G U R E 4 The group differences in functional connectivity in the informative ICs across three groups ( p < .01). For a comparison between X versus Y, silver color: X > Y, cyan color: X < Y.
fewer extreme functional deviations with longer illness durations. As for functional connectivity, the subtype I patients had more abnormalities than healthy controls, and the subtype II patients showed fewer differences compared to healthy controls. Intriguingly, compared with healthy controls, the subtype II patients showed increased and decreased GM. In contrast, the subtype I showed mostly decreased GM. These results herein indicated that the altered brain function of patients became more convergent as the illness duration increased.
However, the range of the illness duration is relatively limited, since the illness of duration of subtype II patients was around 22 months and the duration of subtype I patients was 11 months.
Instead of the traditional case-control studies that tried to localize the differences between two groups, the normative model pro- Ten of the whole brain functional networks are informative at clustering the subtypes of FEDN-MDD. The DMN and frontoparietal brain networks are the most frequently reported abnormal brain network in MDD. A previous study found two neural subtypes of MDD associated with different dysfunctional DMN connectivity patterns (Liang et al., 2020). Another study that used the machine learning approach found that the most discriminative features are distributed in the visual network, somatomotor network, attention network, limbic network, frontoparietal network, and DMN (B. Yan et al., 2020).
Of the two subtypes we identified, when compared to healthy F I G U R E 5 The group differences in gray matter in the informative ICs across three groups (p < .05, cluster size >200). For a comparison between X versus Y, warm color: X > Y, cold color: X < Y. controls, the subtype II patients showed significantly decreased connectivity in DMN. In general, the subtype I patients showed more increased connectivities, compared to healthy controls, especially in the DMN and the limbic systems (the hippocampus-hub, the caudatehub and the thalamus-hub networks, for details, see Figure 4). The increased and decreased functional connectivity in subtype I and subtype II patient could explain the inconsistent results in previous studies. The limbic system comprises key emotional areas that have been reported to be increased- (Jiao et al., 2020), decreased-(Goya-Maldonado et al., 2016, or normally (Dai, Zhou, Xu, & Zuo, 2019) connected in depression. And similar inconsistency was also reported in the within-network FC of the executive control network (ECN), as excessive (Yao et al., 2019) or deficient (Stange et al., 2017) FC have been reported in depression. In a prior study, MDD patients were divided into two biotypes according to the functional connectivity changes in the DMN, one patient subgroup was characterized by hyperconnectivity within the DMN, and one subgroup had decreased connectivity (Liang et al., 2020). In the present study, the pattern was similar but not limited to the DMN, for most brain networks, with one group showing hyper-connectivity while the other showed hypoconnectivity compared to healthy controls. In addition to distinguishing subtype features, some abnormal network patterns (e.g., the decreased functional connectivity between the dorsal and lateral frontal cortex in IC 5) were shared across the two patient subgroups. A previous study compared the functional connectivity of the thalamus in 31 first-episode MDD and 28 healthy controls and found that the patients showed significant abnormal thalamus functional connectivity. This abnormality is correlated with the severity of symptoms (Kang et al., 2018). In the present study, abnormal functional connectivities are mainly found in subtype I patients. In contrast, subtype II patients had a distribution that was less deviated from healthy controls and had relatively longer illness durations. Importantly, the severity of depression and anxiety did not differ between subtypes, indicating that our subtype delineation did not merely separate mild and more severely ill patients. In a recent study, reduced GM volumes were consistently found in two cohorts of clinically diagnosed MDD participants (Hellewell et al., 2019). The structural changes in subtype II patients were more complex than the subtype I patients, since there are both increased and decreased GM abnormalities distributed widely in brain regions in subtype II patients, and most of the abnormalities in subtype I were decreased GM. Several previous studies indicated that the volumes or the FA were reversely correlated with the illness duration (Serra-Blasco et al., 2013;Zhao et al., 2021). In the present study, the subtype II patients have relatively longer illness duration, and more complex brain structural abnormalities. A recently published study divided MDD patients into two groups based on the illness duration, and they found that the network difference in MDD from HCs was confined to those patients with shorter illness durations (Xu et al., 2021). In this study, the difference in illness duration between the two subtypes is not significant enough to get a positive conclusion on the relationship between illness duration and structural/functional brain changes. However, the present results posed an essential opportunity for future studies to investigate structural and functional correlations with illness duration.
Several limitations of this study should be considered. First, our study included treatment-naive patients. The medication's effects on functional connectivity and symptoms did not confound the data analysis. However, it is unclear whether our findings apply to MDD groups with multiple prior episodes and active antidepressant therapies. Second, in certain respects of the method, the ideal choice for the clustering algorithm would be a powerful clustering method capable of detecting clusters of very different types. However, k-means was used in this study, which seeks compact clusters whose time complexity is only linear in the number of data items. Furthermore, the empirical number of ICs should be estimated automatically. Third, this data was from a cross-section study, which lacks treatment information; the value of this clustering on treatment prediction remains to be known.
Forth, the generalizability of the current finding is limited since there are only 14 patients as subtype I in a specific FEDN cohort.
In the present study, we applied the normative model to connectivity network measures and fluctuation coefficients of the brain inherent networks. Combined with the data-driven multivariate clustering method, we recognized the patient cohort into two subtypes.
Two subtypes of patients showed different illness durations and significantly different functional connectivities with GM alterations in several essential brain networks. Our findings provided insight into FEDN depression heterogeneity in biological terms, and represented a promising step toward categorical subtype recognition and personalized treatments.

CONFLICT OF INTEREST
The authors report no disclosures relevant to the manuscript.

DATA AVAILABILITY STATEMENT
The imaging and clinical data used in the present study were obtained from the International Big-Data Center for Depression Research (for details, http://yanlab.psych.ac.cn/IBCDR). The data and code are available upon request from any qualified investigator.

FUNDING STATEMENT
This work was supported in part by Beijing Natural Science Founda-