Major depressive disorder (MDD) is a complex and heterogeneous syndrome with large interindividual variability in symptom profiles, clinical trajectories, and treatment outcomes, affirming the consensus that it is not a singular, unitary disease entity (Goldberg, 2011). The frequent co-occurrence of depression with other psychiatric disorders, such as anxiety, posttraumatic stress disorder, and attention-deficit/hyperactivity disorder, hinders the progress of understanding its biological mechanisms. Traditional analytical methods focus on group averages, operating under the assumption that MDD patients are homogeneous. These methods, typically utilizing tests of statistical significance based on group means, often view interindividual differences as statistical error or noise. This failure to incorporate heterogeneity in statistical models of MDD potentially diminishes the power to detect changes (Buch and Liston, 2021). Therefore, there is an urgent need to conduct studies aimed at understanding the underlying biological mechanisms of heterogeneity in depression, delineating individual-specific mechanisms and treatments.

Early efforts to handle this heterogeneity have attempted to define depression subtypes based on various symptom profiles, and subsequently tested their association with specific neurobiological biomarkers (Maglanoc et al., 2019; Toenders et al., 2020). These categorical definitions have yielded some significant insights into the neurobiology and neural circuitry of depression, particularly in patients who exhibited greater impairment or clinical differences compared to others with the same diagnosis. However, they have also raised several unexpected but important questions. Symptom-based subtypes, for instance, may fluctuate over time. Studies suggested higher stability for psychotic depression and severity-based subtypes but lower stability for melancholic, agitated, and atypical depression subtypes (Coryell et al., 1994; Lamers et al., 2012; Melartin et al., 2004). Moreover, identical symptom profiles may be linked to distinct biological mechanisms, rendering the identification of biological biomarkers for individual-level diagnosis challenging (Beijers et al., 2019).

An increasingly adopted alternative is to dissect heterogeneity using neuroimaging data, which has demonstrated more stable subtyping outcomes than symptom-based approaches (Drysdale et al., 2017). Researchers have employed functional (Liang et al., 2020; Price et al., 2017a, 2017b) and structural (Liang et al., 2019) magnetic resonance imaging (MRI) measures to investigate the heterogeneity and subtypes of MDD. Large samples are often required to model the biological mechanisms of depression heterogeneity (Lynch et al., 2020). MRI offers a key advantage due to its accessibility and the potential for integrating and analyzing data from multiple sites (Schmaal et al., 2020). Moreover, resting-state functional connectivity (RSFC) exhibits high levels of individual specificity and temporal stability (Finn et al., 2015; Gratton et al., 2018), making it a promising tool for investigating heterogeneity within disorders. For example, using RSFC, researchers were able to identify two subtypes in MDD patients that exhibited distinct functional connectivity profiles of the default mode network (DMN) (Liang et al., 2020). One subtype demonstrated increased connectivity while the other subtype exhibited decreased connectivity. These subtypes remained relatively stable across different validation samples (Liang et al., 2020). However, these subtyping studies tend to adopt a categorical approach, overlooking the continuous interindividual variations that exist in the subtypes. Consequently, there is a growing recognition of the need for a dimensional approach that can capture the complex clinical reality of MDD more accurately. To address this need, we employed latent Dirichlet allocation (LDA), a data-driven Bayesian model. This model has been shown to successful in decomposing whole-brain RSFC data to identify dimensional subtypes in autism spectrum disorders (Tang et al., 2020) and to estimate latent atrophy factors based on the gray matter volume (GMV) in Alzheimer's disease (Sun et al., 2019; Zhang et al., 2016) and posterior cortical atrophy (Groot et al., 2020). This approach allows each individual to express multiple disease factors to various degrees rather than assigning them to a single subtype.

In this study, we aimed to provide a detailed characterization of the nature and spatial extent of dysfunctional connectivities in MDD under accounting for significant heterogeneity among individuals with MDD. To reach this goal, we applied the LDA to decompose the RSFC patterns in MDD patients into multiple latent hyper- and hypo-connectivity factors using a large, multi-site sample drawn from the REST-meta-MDD project (Yan et al., 2019). To ensure the robustness of findings, we also tested various analysis strategies (i.e., different preprocessing methods, parcellation schemes, sample sizes, head motion control, education-matched, and illness duration effect estimation). We hypothesized that there were shared and distinct dysfunctional connectivities among factors in MDD and each patient expresses multiple factors to varying degrees.