This study is a secondary analysis of rs-fMRI data of healthy and depressed volunteers recruited starting in 2012 from four medical institutions in Hiroshima, Japan (15). Prior to the administration of any experimental procedure, written informed consent was obtained from all participants. The experiments were carried out in accordance with relevant guidelines and regulations, and all our experimental protocols were approved by the Ethics Committee of Hiroshima University. The mental condition of the volunteers were evaluated using Mini International Neuropsychiatric Interview (M.I.N.I) according to the Diagnostic and Statistical Manual of Mental Disorders, 4th ed. (DSM-IV) criteria (22). Initial screening resulted to 281 healthy participants with no history of any mental or neurological disorder, and 281 diagnosed with major depressive disorder (MDD).

Additional screening of participants was performed to emphasize the boundary between groups in the data set. First, the participants were screened based on their Beck Depression Inventory-II scores (BDI-II; healthy ≤ 13; depressed ≥ 20) (23). Second, participants were removed if they have incomplete data during interview or excessive head movements during fMRI recording. Third, healthy female participants were randomly sampled to match the number of healthy males. Last, the depressed participants were categorized into the melancholic group if they exhibited melancholic features (M.I.N.I, DSM-IV criteria) (22). Otherwise, they are categorized into the non-melancholic group.

In the end, a total of 142 healthy (71 F / 71 M) and 120 depressed (60 F / 60 M) participants were screened and selected for analysis. Within the depressed group, there were 89 melancholic (44 F / 45 M) and 31 non-melancholic (16 F / 15 M) participants. Details on their demographic data are summarized in Table 1.

We recorded resting-state functional Magnetic Resonance Imaging (rs-fMRI) data of the participants using gradient echo planar imaging (EPI) sequences. The imaging device and parameters differ depending on the recording site (15). During recording, participants were instructed to focus on a cross mark displayed on a monitor, and to avoid thinking of anything or falling asleep. The fMRI recording lasted for 5–10 min, depending on the recording site (Table 2).

Data preprocessing for the rs-fMRI data was performed using SPM8 (Wellcome Trust Centre for Neuroimaging, University College London, UK) on Matlab (Mathworks inc., USA). The preprocessing procedure included slice-timing correction, mean image realignment, normalization and resampling through segmentation of structural image aligned with the mean functional image, and smoothing with an isotropic 6 mm full-width half-maximum Gaussian kernel (24). The potential confounding effects (i.e., the temporal fluctuations of the white matter, cerebrospinal fluid, and the global signal, as well as six head motion parameters) were linearly regressed out from the fMRI time series to remove the physiological noise and motion artifacts (25, 26). These components were determined by the T1 images, which were simultaneously recorded with the rs-fMRI data. Finally, head motion artifacts were scrubbed from the functional images based on the relative changes (i.e., translational displacements along X, Y, and Z axes, and rotational displacements of pitch, yaw, and roll) between the image frames through time, with a frame-wise displacement (FD) threshold = 0.5 mm (27).

When selecting the regions of interest (ROIs) for analysis, we considered functional brain networks that were associated with depression. Time series data were extracted from 74 ROIs according to the 90-ROI Shirer Brain Atlas (28). Some ROIs were excluded due to either unreliable images during recording (e.g., in cerebellum area) (29), or insufficient number of ROIs of corresponding functional network (e.g., primary and higher visual networks). When extracting the time series data, the global mean and confounding effects of both cerebrospinal fluid (CSF) and white matter were linearly regressed out as part of the preprocessing step. A total of 12 distinct functional brain networks were analyzed in this study. For the brevity of this paper, we highlighted most of the results for three main networks: basal ganglia network (BGN), dorsal default mode network (DDMN), and left executive control network (LECN). These networks have been associated with melancholic depression (14, 30), and melancholic symptoms such as anhedonia (31) and rumination (32). The list of anatomical ROIs for all of the analyzed networks networks are listed in Table 3 and Supplementary Table 1.

To compare the brain dynamics of melancholic and non-melancholic participants, we implemented the Energy Landscape Analysis (ELA) method which utilizes the Pairwise Maximum Entropy Model (P-MEM) (Figure 1) (36). In this study, ELA was conducted separately for each of the 12 distinct functional networks (see Supplementary Table 1).

In ELA, we began by combining the rs-fMRI data of all participants in a group (healthy, non-melancholic, melancholic) to create group-level dynamics model. In each group, the concatenated time series data were extracted from the corresponding ROIs of the chosen functional network (Figures 1A,B). Then, the time series data were converted to binarized signals, using the average signal value as the threshold (Figure 1C). At any time point, binarized signals were either 1 (active) or 0 (inactive). The threshold value was computed for each ROI across the combined data of the group participants.

Binarizing the signals allowed us to encode the time signals into brain state sequences, wherein a brain state was defined by the activity pattern, i.e., active and inactive ROIs, at a given time (Figures 1C,D). Given brain network with n ROIs, there are 2ⁿ possible states. For each brain state σ, we computed the empirical probability p(σ),

where n_σ is the number of occurrences that the state σ appeared in the time series, and T is the total length of the time series (Figure 1D).

From the brain state probabilities, we defined the state energy using Boltzmann distribution (37),

where the sum in denominator is taken over all 2ⁿ possible states (each state is denoted by σ′). In P-MEM, the state energy E(σ) in Equation 2 was restricted in the quadratic form, and expressed as:

with h=[hi]i=1n and J=[Ji,j]i,j=1n corresponding to the individual ROI activity and the pairwise ROI interaction, respectively. As implied by these equations, more active ROIs correspond to higher h and J, which leads to lower (more negative) energy (Equation 3), and higher occurrence probability (Equation 2) (36).

We then built the P-MEM from the model parameters h & J by maximizing their likelihood, as defined by

for the entire time series of length T.

With the P-MEMs for healthy, non-melancholic and melancholic groups, we analyzed the energy landscapes and transition dynamics of heterogeneous depressed brain networks (Figure 1E). The ELA features analyzed in this study are summarized in Supplementary Table 2.

To start with, basin state dendrograms were constructed by finding the basin states and their clusters. Basin states (or basins) are states with the lowest energy relative to their neighboring states. States are said to be neighbors if they differ in only one active/inactive region. If neighboring states have lower energy barriers relative to other states, they become part of the closest basin cluster (38). Basins are the core of energy landscapes since these states are presumed to be the most stable. Djikstra's algorithm was performed to search for the basins, and construct the leaves (basins) and branches (energy barriers) of the dendrogram. In Figure 1E: Energy Barriers, the energy barrier between states a and b are highlighted (red: a → b; blue: b → a). Here, a → b has higher barrier, and thus has lower probability of occurring (38).

Energy landscapes were then constructed to depict the energy level and cluster size of the basins. In Figure 1E: Energy Landscape Analysis, a 3D schematic diagram depicts the basins in an arbitrary state space. Here, basin clusters are plotted as concentric circles, where the basin is at the center, and neighboring states are on circles with radius equivalent to the their distance to the basin (i.e., a state that differs in one region with the basin has distance of 1, state the differs in two regions has distance of 2, etc.). Thus, each circle represented the states that were equidistant to the basin. Then for each circle (including the center), its depth was equivalent to the energy of the state with lowest energy (Equation 3) in that circle.

Since energy landscapes modeled only the group-level brain dynamics, we opted to analyze the brain dynamics of individual participants by computing the occurrence frequency of basins f(A) on individual fMRI time series. The occurrence frequencies of all basins (major and minor) are then summed for each group and for each network.

We also analyzed the transition dynamics on the energy landscapes. For each network, we selected two major basins as the basins with the lowest energy, A₁ and A₂; with P₁ as clustered states for A₁, P₂ for A₂. Then we computed the following participant-level transition rates,

where n(U → V) is the number of times the participant's time series entered state U and arrived at state V; TR(A → A′) is the major transition rate for major basins; TR(P → P′) is the peripheral transition rate for clusters of major basins. Moreover, we computed the staying rates of the states,

where this time we counted the number of times the participant's brain activity stayed on the major basins (SR(A → A′)), or the periphery (SR(P → P′)). The transition rates were used to define the Traveling Score, which measured the amount of times the brain successfully traveled from one major basin to another (Equation 10, Figure 5A). Similarly, the staying rates were used to define the Lingering Score, which measured the amount of times the brain lingered on a basin or along its periphery (Equation 11, Figure 5A),

Note that the traveling score (Equation 10) is also known as efficiency score, which is the index of ease of transitions (38).

One-way Analysis of Variance (ANOVA) was performed to find significant group differences (healthy vs. non-melancholic vs. melancholic) in participants' age, IQ, BDI-II and SHAPS (Snaith-Hamilton Pleasure Scale) scores. Bonferroni correction was applied to compensate for multiple group comparisons. Similarly, Chi-squared test was performed for group differences in sex. Group differences are presumed to be statistically significant at p < 0.05 for both types of statistical analyses. One-way ANOVA with Bonferroni correction was also applied to test for statistical significance of results when comparing individual- and group-level differences in energy landscape characteristics such as basin frequency and transition rates. These results were also verified by Kruskal-Wallis test to cope with non-normally distributed data as well.

Finally, we investigated the correlation between depressive symptoms and basin characteristics. Two criteria were used for analyzing depressive symptoms: BDI-II for depression diagnosis (23), and SHAPS for assessing anhedonia (39). Multivariate linear regression model were fitted on the data, with p < 0.05 deemed as statistically significant relationship between the variables.

We found significant differences between healthy, non-melancholic, and melancholic groups in terms of BDI-II and SHAPS score (Table 1). For both criteria, the melancholic group had higher average scores (p < 0.05, Healthy vs. Melancholic; p < 0.001, Non−melancholic vs. Melancholic), followed by non-melancholic group (p < 0.001, H vs. N). The results were consistent with many previous reports (40–42).

In terms of other demographical factors, we found no significant group differences in sex distribution, age, and IQ levels of the participants. However, there was significant difference between the number of melancholic and non-melancholic participants recruited on four different sites (Table 1).

The basin dendrograms of each group for all network revealed differences in major basins and energy levels. The dendrograms for Basal Ganglia Network (BGN), Dorsal DMN (DDMN), and Left ECN (LECN) are shown in Figure 2. Here, the major basins (purple) were deeper (i.e., had relatively lower energy) than the minor basins (green).

In both BGN and LECN, the non-melancholic and melancholic groups had deeper major basins than the healthy group. Remarkably, a different pattern can be observed in the DDMN, where the non-melancholic group had shallower major basins than the healthy and melancholic groups.

Furthermore, some major basins were unique to specific groups. For example, non-melancholic group had a unique major basin in DDMN, where the left and right angular gyri (LAG, RAG), posterior cingulate cortex (PCC), and precuneus regions were activated. Similarly on LECN, the non-melancholic group also had a unique major basin, where left inferior frontal gyrus (LIFG), orbitofrontal gyrus (OFG), left inferior temporal gyrus (LITG), middle temporal gyrus (MTG), and left thalamus regions were activated. Complete results for basin energies are reported in Supplementary Table 3.

Lastly, the major basin pairs had antisynchronized activations (Figure 2). In BGN, the Left and Right Inferior Frontal Gyrus (LIFG, RIFG) were antisynchronized with the Left and Right Thalamus, Caudate, Putamen and Pons. In DDMN, the Left and Right Angular Gyrus (LAG, RAG), Right Superior Frontal Gyrus (RSFG), Posterior Cingulate Cortex (PCC) and Precuneus were synchronized with each other. A similar pattern can be observed in LECN, where the Left Inferior Temporal Gyrus (LITG), Middle Temporal Gyrus (MTG) and Left Thalamus were synchronized.

Energy landscapes are then constructed based on the basin information. In Figure 3, different energy landscapes can be observed from each group in DDMN. The melancholic group has the most number of basins (n = 6), followed by healthy group (n = 5), then non-melancholic group (n = 4). Because of this shortage of basins, the non-melancholic has relatively larger major basins (s_A = 98.6%), compared to healthy (s_A = 92.0%) and melancholic groups (s_A = 89.4%). Note that, as mentioned previously, the non-melancholic group has a unique major basin (A₂). Also, there is a hill on basin cluster B₂ of the non-melancholic group, which may be inferred as highly unstable cluster of states. This hill state cluster is absent in the other groups. We compared the basin sizes (number of states clustered to the basin) of each group for each network but found only minimal differences (e.g., p = 0.08, N vs. M, DDMN). Complete results for basin sizes are reported in Supplementary Figure 1 and Supplementary Table 3.

Significant group differences were found on the basin frequencies across all selected networks. The most notable results are summarized in Figure 4. In some networks, there was a natural decreasing trend from healthy to melancholic group (e.g., p < 0.005; H vs. M, N vs. M; Auditory Network). In other networks, an opposite, increasing trend appeared (e.g., p < 0.005; H vs. N, H vs. M, N vs. M; Anterior Salience Network, DDMN). And for the rest, the non-melancholic group had either the lowest (e.g., p < 0.005; H vs. N, H vs. M, N vs. M; BGN, LECN) or highest frequency (e.g., p < 0.005; H vs. N, H vs. M, N vs. M; Posterior Salience Network). Similar results for the other functional networks are summarized in Supplementary Table 4.

Depression heterogeneity was also evident when comparing the transition dynamics in some networks. As shown in Figure 5B, there were significant differences between non-melancholic and melancholic groups in traveling scores (p < 0.05; N vs. M; DDMN), and lingering scores (p < 0.01; N vs. M; DDMN, LECN) of individual participants. There was also significant group differences in the lingering scores between healthy and non-melancholic (p < 0.01; H vs. N; DDMN), and healthy and melancholic groups (p < 0.01; H vs. M; DDMN, LECN). Furthermore, there is significant increase in lingering scores of melancholic group in LECN (Lingering = 0.61 ± 0.06) as compared to non-melancholic (Lingering = 0.55 ± 0.06; p < 0.005) and healthy groups (Lingering = 0.53 ± 0.06; p < 0.005). Results for other functional networks are summarized in Supplementary Table 5.

Our first goal was to show heterogeneity in our depressed participants. Depression heterogeneity has plagued research advancement on the diagnosis and treatment of the disorder (47–49). Considerable efforts has been made toward drawing boundaries within the currently accepted definition of major depressive disorder (47, 48), suggesting that investigations on individual depressive symptoms and their interactions may lead to new discoveries (49).

Melancholia is recognized as a specifier for MDD, characterized by symptoms that overlap with MDD such as anhedonia, excessive guilt, psychomotor disturbance, cognitive impairment, weight loss, and worse mood in the morning (2, 50). Although melancholic depression is not anymore considered as a depression subtype by DSM-5 (which instead uses the technical term depression with melancholic features) (50, 51), its clinical characteristics distinguishes it from non-melancholic and atypical depression (43). Our analysis results confirm such distinctions, as observed from significant differences between melancholic and non-melancholic groups in BDI-II and anhedonia scores (Table 1); basin frequencies on 12 networks (Supplementary Table 4); and transition dynamics scores in DDMN, LECN, PSN, RECN, VDMN, and VSPN (Supplementary Table 5).

Next, we sought meaningful relationship between energy landscape signatures and depression severity. Melancholic depression is associated with greater symptom severity in comparison to depression with non-melancholic features (40, 43). This association was evident in our participant-level analysis of BDI-II scores, where the melancholic group scored significantly higher than the other groups (Table 1).

Readers may wonder if energy landscape features can be a proxy to evaluate some symptoms associated with melancholic features. To check this possibility, we performed regression analysis on basin frequencies and depression severity score (BDI-II). But despite having clear boundaries among the groups, no significant correlation were found in any of the networks (Supplementary Figure 2A).

The order of increasing severity from healthy, to non-melancholic, to melancholic depression may seem intuitive. However, some studies argue against the severity-based categorization of melancholic depression (52, 53), and do not identify melancholia as a “more severe” form of depression (52). This is evident in Supplementary Figure 2A, where the BDI-II scores for non-melancholic and melancholic groups coincide. In fact, we can even argue that the basin frequency of our energy landscape model were better than the depressive scores in discriminating melancholic depression (Figure 4).

Following depression severity, we then tested if the correlation with basin frequency would also hold true for anhedonia. Anhedonia is one of the overlapping symptoms for both melancholia and depression (2, 42). This effect was observed from the statistical analysis of the subject demographics, where melancholic group has significantly higher SHAPS score for anhedonia (Table 1).

Similar to BDI-II scores, there were no significant correlation between SHAPS scores and basin frequencies (Supplementary Figure 2B). However, a closer inspection of results in the DDMN reveals a plausible connection between these features. Studies report that the decrease in functional connectivity of dorsomedial prefrontal cortex (dmPFC) with posterior cingulate cortex / precuneus (PCC/PCUN) is related to depression severity and anhedonia (33, 44).

Looking at the major basins found by our model in DDMN (Figure 2), the major basin pairs for all groups have antisynchronized activations of mPFC and PCUN (S₂₇₉:{mPFC = 1, PCC/PCUN = 0}, H, N, M; S₂₃₂:{mPFC = 0, PCC/PCUN = 1}, H, M; S₁₆₈:{mPFC = 1, PCC/PCUN = 0}, N). This decoupling between mPFC and PCC/PCUN may be attributed to impaired reward anticipation, which is a key feature of anhedonia (44). Although current research has not yet established the direct connection between dynamic functional connectivity and energy landscape model (38), our results provide insight to the potential connection between these models.

Lastly, we investigated the possible presence of ruminative states in the energy landscapes. Rumination is the tendency to dwell on the same, usually negative, thoughts for prolonged periods of time. Currently, rumination is not recognized a diagnostic feature for melancholic depression. However, since past studies have noted the significant correlation of rumination to both melancholy (45) and depression (54), we explored the energy landscape characteristics to provide more evidences on this correlations. Ruminative thinking is marked by increased connectivity in DMN (55), and decreased activity in left dorsolateral prefrontal cortex (left dlPFC; equivalent to LECN) (46). In terms of our energy landscape model, we define dynamic activity as being able to transition from one basin to another. This is in contrast to static activity, which tends to stay at the same basin. We derived two measures for brain activity: traveling score (Equation 10) for dynamic activity, and lingering score (Equation 11) for static activity. Thus, high traveling score and low lingering score would imply a more active network.

The significant increase of lingering scores of the melancholic group in LECN (Figure 5) could imply that the melancholic group has greater tendency to be “stuck” in a basin, or within its cluster. This may be related to the decreased activity in left dlPFC found in rumination (46). On the contrary, in the dorsal part of DMN (DDMN), the lingering scores of melancholic group (Lingering = 0.55 ± 0.07) are significantly lower than healthy (Lingering = 0.58 ± 0.06; p < 0.01) and non-melancholic groups (Lingering = 0.63 ± 0.05; p < 0.005). Although this may support previous studies showing increased connectivity in DMN during ruminative thinking, analyzing the lingering scores in ventral part of DMN (VDMN) leads to contradicting results (55). Nevertheless, the traveling scores in DDMN are significantly higher in melancholic group (Traveling = 0.10 ± 0.09) than in non-melancholic group (Traveling = 0.05 ± 0.06; p < 0.05). For LECN and VDMN, there were no significant differences in traveling scores between groups.

Unlike typical functional connectivity (FC) models that are based on correlation between two regions (i.e., using Pearson's correlation coefficient), the ELA model does not assume that pairwise interactions of regions are independent from each other (17). This allows our model to more accurately capture the global activity patterns that may possibly be overlooked by FC-based models (17). Some FC-based models try to work around this limitation by using a different FC metric (such as precision matrix), or by introducing a sliding window in the analysis (also known as dynamic FCA) (56).

Hidden Markov Model (HMM) is another common brain dynamics model for resting-state fMRI. In comparison to ELA which utilizes pairwise MEM, HMM is more complex, and may be more expressive (57). Thus, in the future, our ELA results can serve as a baseline for deeper analysis on brain dynamics of melancholic depression using models such as HMM.

Furthermore, our present study have other limitations. ELA assumes a memoryless process, such that predicting the next state only depends on the current state, and not the past states (38). This saves us memory space and processing time (58), but this only works on systems at equilibrium, which typically requires huge amount of data (59). To address this problem, first, we decided to use resting-state fMRI data since these are expected as equilibrium states (60). Then, we combined the fMRI time signals of all participants in a group, to essentially produce “long” temporal data. Since there are substantially more healthy participants, it is possible that the model for the healthy group has converged closer to equilibrium than the two depressed subgroups.

Although concatenation of individual fMRI data should only have minimal effect to our analysis due to the memorylessness aspect of ELA (61), it should be noted that combining participants fMRI data could still potentially introduce bias due to BOLD signal differences across participants. We addressed this by using individual subjects' mean fMRI BOLD signals as threshold for binarizing the signals. However, it is recommended to further investigate the implications of data concatenation, and the ways to mitigate the intensity differences in individual fMRI.

Our selection of ROIs relies on the parcellation method [Shirer Atlas, (28)] and availability of data. As such, we removed some ROIs due to unreliable or insufficient data. Inclusion of these ROIs might affect the results. Choosing a different parcellation method might also produce different results. Thus, testing the robustness and consistency of our model on complete sets of ROIs or with different parcellation methods would be a reasonable step in the future.

It is often pointed out that the site difference could be a potential confounding factor for the observed BOLD signals (62). In fact, when we applied two-way ANOVA considering two factors of participant groups and site ID to basin frequencies, all main effects and their interaction were statistically significant (Supplementary Table 4). While this still supports that basin frequencies across different networks could differentiate participants suffering from melancholic and non-melancholic MDD in a population level, its reliability as a biomarker for the personalized medicine should never be overestimated. Thus, the potential site bias is a limitation in the present study and we should address the issue in our future work.

Even though we hypothesized a priori that the energy landscape of melancholic depression will be different from non-melancholic due to depression heterogeneity, our analyses and interpretations of the results were done post-hoc. For statistical analyses, we applied Bonferroni correction to compensate for multiple comparison tests that would increase the risk of Type 1 error (63).

Our discussion on the relation between melancholia and rumination is more suggestive than conclusive since we lack rumination quantifiers to analyze with our models. In the future, it would be crucial to record ruminative tendencies of participants using standardized measures such as Ruminative Response Scale (RRS) (64).

Finally, we would like to emphasize that melancholia, in itself, is also heterogeneous (50). It is characterized by multiple symptoms, and thus its severity and symptoms may vary from one person to another. Current biomarker models focus on the distinguishing features of melancholic depression on a group level (9, 10, 14). ELA allows us to study not only the group-wide brain dynamics of melancholic depression, but also the individual nuances in brain states and functional region interactions. In the future, it would be beneficial to study more of these individual differences so we can directly confront the issue of heterogeneity.

Melancholic depression is a debilitating disorder that robs a person the pleasure and excitement from activities they used to enjoy. In this study, we developed an energy landscape model to better understand the brain dynamics involved in melancholic depression. Relative to healthy and non-melancholic groups, the melancholic group showed significant differences on basin energy, basin frequency, and transition dynamics in several functional networks. Moreover, possible connections were traced between major basins and depressive symptom scores such as depression severity and anhedonia. Taken together, these results suggest that melancholic depression is a distinct disorder, and thus should be diagnosed and treated with utmost precision and care.