Manifold alteration between major depressive disorder and healthy control subjects using dynamic mode decomposition in resting-state fMRI data

Background The World Health Organization has reported that approximately 300 million individuals suffer from the mood disorder known as MDD. Non-invasive measurement techniques have been utilized to reveal the mechanism of MDD, with rsfMRI being the predominant method. The previous functional connectivity and energy landscape studies have shown the difference in the coactivation patterns between MDD and HCs. However, these studies did not consider oscillatory temporal dynamics. Methods In this study, the dynamic mode decomposition, a method to compute a set of coherent spatial patterns associated with the oscillation frequency and temporal decay rate, was employed to investigate the alteration of the occurrence of dynamic modes between MDD and HCs. Specifically, The BOLD signals of each subject were transformed into dynamic modes representing coherent spatial patterns and discrete-time eigenvalues to capture temporal variations using dynamic mode decomposition. All the dynamic modes were disentangled into a two-dimensional manifold using t-SNE. Density estimation and density ratio estimation were applied to the two-dimensional manifolds after the two-dimensional manifold was split based on HCs and MDD. Results The dynamic modes that uniquely emerged in the MDD were not observed. Instead, we have found some dynamic modes that have shown increased or reduced occurrence in MDD compared with HCs. The reduced dynamic modes were associated with the visual and saliency networks while the increased dynamic modes were associated with the default mode and sensory-motor networks. Conclusion To the best of our knowledge, this study showed initial evidence of the alteration of occurrence of the dynamic modes between MDD and HCs. To deepen understanding of how the alteration of the dynamic modes emerges from the structure, it is vital to investigate the relationship between the dynamic modes, cortical thickness, and surface areas.


Introduction
The World Health Organization has reported that approximately 300 million individuals suffer from the mood disorder known as major depressive disorder (MDD).MDD gives rise to psychological symptoms, such as despondent moods and negative cognitions, as well as physical symptoms, such as sleep disturbances and fatigue in mild cases, and even suicide in severe cases (1).Neurotransmitter reuptake inhibitors, such as selective serotonin reuptake inhibitors and transcranial magnetic stimulation through electrical stimulation, have been employed in the treatment of MDD (2-4).Although these treatments are effective, there are patients whose depressive symptoms improve only partially or not at all (5).Therefore, the mechanisms underlying MDD need to be elucidated.
Non-invasive measurement techniques have been utilized to reveal the mechanism of MDD, with resting-state functional magnetic resonance imaging (rsfMRI) being the predominant method (6).To evaluate dynamic changes in blood oxygenation level-dependent (BOLD) signals using rsfMRI, static functional connectivity (sFC), dynamic functional connectivity (dFC), and energy landscape (EL) were employed as indices to portray the dynamics of whole-brain networks.sFC captures the static relationships of spontaneous fluctuations that represent correlations over the entire duration (7,8), whereas dFC captures time-resolved spontaneous fluctuations in which functional connectivity (FC) changes over a short time (9)(10)(11).Evaluation of the static and dynamic relationships of spontaneous fluctuations in the whole-brain network has revealed that MDD exhibits abnormal connections in FC, such as the default mode network (DMN), control executive network (CEN), and salience network (SN) when compared with healthy controls (HCs) (12)(13)(14)(15)(16). Analyzing sFC involves calculating the correlation between two independent regions for all pairs (17).Even if a pair of regions is not directly structurally interconnected, their sFC can exhibit a strong correlation if both regions receive input from a third region (18).Hence, it is imperative to simultaneously represent the dynamics of whole-brain networks based on neural activity across multiple regions.This is where EL emerges, which utilizes a pairwise maximum entropy model to represent the dynamics of the whole-brain network in terms of the activity within each region and the interactions between two or more regions (19).Moreover, by defining the functional network between subjects in terms of energy, it is possible to evaluate the transition from one stable state to another through the unstable states.Notably, MDD tends to sink to specific states, and it is difficult to transition from one stable state to another compared to HCs (20).Although EL excels in stability analysis across subjects, some issues require prior assignment of a functional network to each region and binarization of BOLD signals.In common with sFC, dFC, and EL, analyzing components of the BOLD signal above 0.1 Hz is a challenging problem.Therefore, in terms of interactions across multiple regions, a methodology is required to evaluate the sinking into specific states under conditions free from functional network assignment and binarization.
The dynamic mode decomposition (DMD) is a data-driven and equation-independent approach for analyzing fluid dynamics (21).DMD calculates eigenvectors and corresponding eigenvalues of the approximate linear transformation expressing the time evolution of multidimensional time-series data.Eigenvectors were called dynamic modes (DMs) representing coherent spatial patterns and the corresponding eigenvalues were called discrete-time eigenvalues representing the frequency and time evolution such as growth and decay.In other words, multiple coherent DMs coexist at a certain time in multidimensional time-series data and corresponding temporal characteristics are identified.EL analysis assigns a functional network to each region, binarizes the BOLD signal, fits it with a Boltzmann distribution, determines relationships between activity patterns and energy, and assigns one state on EL at a certain time in multidimensional time-series data (22).Here, since the BOLD signals exhibit wave superposition, it is necessary to analyze stability under conditions where multiple states coexist at a certain time.DMD was successful and recent studies have applied DMD to BOLD signals, a type of fluid that exhibits nonlinear spatiotemporal changes (23)(24)(25)(26).This study applied DMD to the BOLD signals across all frequency bands of HCs and MDD.Subsequently, the spatial patterns, frequencies, and temporal changes across all subjects were analyzed in terms of stability.
Analysis of a large dataset of psychiatric disorders based on rsfMRI (27) using DMD revealed that the number of DMs associated with MDD decreased in visual networks (VN) and SN, while it increased in DMN and sensory-motor networks (SMN) when compared to HCs.Interestingly, DMs' differences between MDD and HCs were identified not only within the 0.01-0.1 Hz range in standard rsfMRI analysis but also extending beyond 0.1 Hz.Applying t-distribution stochastic neighbor embedding (t-SNE) (28) to DMs enables the disentangling of the intricate curved surfaces spanned by DMs into a two-dimensional manifold, allowing for the evaluation of stability across subjects.Subsequently, DMs resembling resting-state networks (RSNs) were identified by evaluating the probability density ratio between HCs and MDD using a two-dimensional manifold.The amplitudes of the DMs resembling the VN and SN were similar to the spatial patterns associated with cortical thickness and surface area abnormalities in MDD (29).

Materials and methods
In this study, we applied DMD to BOLD signals and devised a method for extracting DMs based on the probability density ratio between HCs and MDD on two-dimensional manifolds using t-SNE (Figure 1).First, the BOLD signals of each subject were transformed into DMs representing coherent spatial patterns and discrete-time eigenvalues to capture temporal variations using DMD.Second, all the DMs were disentangled into a twodimensional manifold using t-SNE.Finally, density estimation and density ratio estimation were applied to the two-dimensional manifolds after the two-dimensional manifold was split based on the HCs and MDD.The results revealed that MDD tended to sink into specific DMs in contrast to HCs.

Dataset
We used the Japanese Strategic Research Program for the Promotion of Brain Science (SRPBS) dataset (27) (https:// bicr.atr.jp/decnefpro/data/),along with additional datasets obtained from various projects.Supplementary Table 1 describes the protocols at each site, and Supplementary Table 2 describes the subject information at each site.
The datasets were collected from the Center of Innovation at Hiroshima University (COI) and the University of Tokyo (UTO), Hiroshima Kajikawa Hospital (HKH), Hiroshima Rehabilitation Center (HRC), Hiroshima University Hospital (HUH), and Yamaguchi University (UYA).COI and UTO follow the unified protocol but HKH, HRC, HUH, and UYA follow non-unified protocols.The total number of HCs and MDD was 543 and 302, respectively, with Beck Depression Inventory-II (BDI-II) scores of 7.5 ± 6.3 and 28.1 ± 10.5, respectively.

BOLD signals preprocessing
BOLD signals were preprocessed using fMRIPrep version 1.0.8(http://fmriprep.readthedocs.io/en/1.0.8/workflows.html)(30).The first 10 s of the data were discarded to allow for T1 equilibration.The preprocessing steps included slice-timing correction, realignment, coregistration, distortion correction using a field map, segmentation of T1-weighted structural images, normalization to Montreal Neurological Institute space, and spatial smoothing with an isotropic Gaussian kernel of 6 mm full width at half maximum."Fieldmap-less" distortion correction was performed for the test dataset due to the lack of field map data.

Preprocess of ROI time series for DMD
It is necessary to mitigate the effects of the protocols and physiological noise.BOLD signal extraction was performed using Glasser's 360 regions of interest (ROI) (31), which excluded the cerebellum and contained little white matter.Overview of the analysis procedure.First, each subject's blood oxygenation level-dependent (BOLD) signals were extracted using Glasser's 360 regions of interest (ROI).Second, the BOLD signals were decomposed into dynamic modes (DMs) and discrete-time eigenvalues using the onestacked time-delay coordinates dynamic mode decomposition (tdcDMD).Third, all DMs were disentangled into the two-dimensional manifold using t-distributed stochastic neighbor embedding (t-SNE).Fourth, density estimation was performed to visualize the features that major depressive disorder (MDD) sink into the specific DMs compared to healthy controls (HCs).Finally, density ratio distributions between HCs and MDD were calculated using relative unconstrained least-squares importance fitting (RuLSIF).Detrending was applied to eliminate long-term variations, and BOLD signals were normalized using z-scores to mitigate the effects of the protocols.When analyzed using the DMD, the frequencies were computed for each DM.Therefore, band-pass filtering was not applied.
Confounding factors must be removed when extracting BOLD signals.The fit _transform function was applied to remove confounding factors for the 12 regression parameters (six motion parameters, average signals over the whole brain, and five anatomical CompCor components).

One-stacked time-delay coordinates DMD
BOLD signals were decomposed into DMs and discrete-time eigenvalues.Time-delay coordinates DMD (tdcDMD) is a method used for decomposing standing waves into spatiotemporal patterns with high accuracy (21); tdcDMD was performed using the dmd.py function in the DMD toolbox (https://github.com/erichson/DMDpack).As described in a previous study (26), the BOLD signals of each subject were converted into DMs.As shown in Equation 1, the BOLD signal matrix X was composed of rows representing the number of ROI, N roi and columns representing the number of measurements, N T .
where x k ( ∈ R Nroi ) represents the BOLD signals at time k.The following matrices were constructed from the BOLD signal matrix X as shown in Equations 2, 3.
where X 2 represents the matrix with X 1 shifted back one observation.Subsequently, x k+1 was stacked on x k as shown in Equations 4, 5.
X 2aug was predicted using X 1aug so X 2aug ≈ AX 1aug .
where the dagger represents the generalized inverse.Singular value decomposition was applied to X 1aug .
where U, S, and V represent the left singular, singular value, and right singular matrices of X 1aug , respectively.As shown in Equation 8, the matrix A is rewritten by substituting Equation 7 into Equation 6.
the proper orthogonal decomposition was applied to A.
then eigen decomposition was applied to Ã.
where W and L represent the eigenvector and eigenvalue matrices of Ã, respectively.X 2aug VS −1 was multiplied from the left in Equation 10 and Equation 9 was substituted into Equation 10.
where Ã is the similar matrix of A, so they have the same eigenvalue matrix L but different eigenvector matrices.In comparing Equations 10, 11, X 2aug VS −1 W can be regarded as the eigenvector matrix of A. Finally, the eigen decomposition of A was reconstructed using W and L and the dynamic mode matrix F was calculated as shown in Equation 12.
the i-th column of F, which we denote by f i ( ∈ C 2N roi ), is the ith eigenvector of A. The i-th diagonal element of L, which we denote by l i ( ∈ C), is the i-th eigenvalue of A. The phase and amplitude of l i mean the frequency and decay rate of the corresponding mode.The frequency f i corresponding to the dynamic mode f i and the eigenvalue l i is described as following Equation 13.
where Dt, ln( • ) and imag( • ) represent the temporal resolution in each protocol, natural logarithm, and the imaginary part of a complex number.

Two-dimensional manifold with t-SNE
When analyzed using the DMD, pairs of DMs with identical amplitudes but antiphases emerged.Moreover, DMs representing brain states describe intricate curved surfaces in a multidimensional space.In a previous study (26), the modified K-means clustering algorithm was applied to DMs and treated DMs with identical amplitudes and antiphases.However, this approach failed to disentangle intricate curved surfaces in a multidimensional space.Hence, this study employed t-SNE (28) to disentangle the intricate curved surfaces spanned by DMs.
The initial 360 rows, which are inherently independent of the 720 rows of the DMs, were used to employ a one-stacked tdcDMD.Subsequently, the DMs were separated into their real and imaginary components, stacked together, and applied to t-SNE.When t-SNE was applied to all DMs of both HCs and MDD, the perplexity varied from 30 to 10,000.A value of 2,000 was visually selected to achieve maximum separation between peaks within the two-dimensional m a n i f o l d w h i l e k e e p i n g r a n d o m _ s t a t e fi x e d .T h e sklearn.manifold.TSNE function in Python was employed, with all parameters set to their default values except perplexity and random_state.

Kernel density estimation
It is crucial to select the optimal perplexity at which the peaks within the two-dimensional manifold achieve maximum separation.Hence, we separated the peaks by performing a kernel density estimation on a two-dimensional manifold.The formula for estimating the probability density r at a given point y, estimated from points x i (i = 1, 2, …, n) of DMs on the two-dimensional manifold is expressed as following Equation 14: where kernel K is the Gaussian kernel and bandwidth h is set to the Scotts factor.Scipy.stats.gaussian_kdefunction in Python was used (32).

Kernel density ratio estimation
The probability density was estimated using kernel density estimation on the two-dimensional manifolds obtained by applying t-SNE.Consequently, the distinction between HCs and MDD was revealed as a different balance in the proportion of DMs rather than the emergence of unknown DMs.Hence, we estimated the probability density ratio between HCs and MDD using a relatively unconstrained least-squares importance fitting (RuLSIF) (33).In terms of estimation accuracy, it is more precise to directly estimate the density ratio between HCs and MDD than to indirectly estimate the density ratio by estimating HCs and MDD's densities separately and dividing HCs and MDD's densities.To improve the estimation accuracy, various methods have been developed to directly estimate the density ratio without going through the density estimation process.RuLSIF was chosen for this study because its Python code is publicly available and its calculation speed is fast.
To estimate the density ratio of the area where the HCs' density was higher than the MDD's density, the HCs' manifold was used as the denominator, and the MDD's manifold was used as the numerator.To estimate the density ratio for the area where the MDD's density was higher than the HCs' density, the MDD's manifold was used as the denominator, and the HCs' manifold was used as the numerator.
2.8 Plotting dynamic modes, histogram of frequency, and discrete-time eigenvalues greater than 95% significance level Kernel density ratio estimation was used to calculate the probability density ratio between HCs and MDD.However, the specific regions exhibiting significant differences in terms of density ratio between HCs and MDD remain unknown.To solve this problem, permutation tests were performed to clarify areas higher than the 95% significance level and to plot the mean amplitude and phase of the DMs, a histogram of frequency, and discrete-time eigenvalues within the significant areas.
First, we randomized the labels of the HCs and MDD in a twodimensional manifold.Second, with fixed parameters (a, s , h) = (0, 1:0, 0:01), RuLSIF was performed to calculate the maximum peak value, repeating this process 100 times.Third, we applied the density-based spatial clustering of applications with noise (DBSCAN) (34) to cluster points within areas that exhibited maximum peak values higher than the 95th percentile.Finally, we plotted the mean amplitudes and phases of the DMs, frequency histograms, and discrete-time eigenvalues l associated with each cluster.For the density ratios p MDD (x)=p HCs (x) and p HCs (x)=p MDD (x), the DBSCAN parameters were set as (eps, min samples) = (1, 100) and (0:15, 300), respectively.Points that were not assigned to a cluster were excluded.

Applying t-SNE, density estimation, and density ratio estimation to the DMs
First, the two-dimensional manifold was calculated by applying t-SNE to all DMs across all subjects and was visualized after separating the HCs and MDD (Figure 2A: HCs, B: MDD).Second, the perplexity was varied from 30 to 10,000 and consequently set to 2,000 to maximally separate the peaks in the two-dimensional manifold.Finally, kernel density estimation was performed to clarify the distribution features exhibited by the twodimensional manifold (Figure 2C: HCs, D: MDD).
In the HCs, the peaks displayed a relatively uniform distribution (Figure 2C).Conversely, in the MDD group, the peaks exhibited a bias toward the upper right, lower left, and central areas (Figure 2D).In other words, MDD tended to sink more into specific DMs than HCs.In addition, the edge of the MDD manifold appeared slightly wider than that of the HCs manifold at the elliptical periphery.To assess these features, density ratio estimation was performed by applying RuLSIF to the twodimensional manifolds.

DM's features in the clusters
The density ratio was calculated using HCs as the denominator and MDD as the numerator (Figure 2E).Similarly, the density ratio was calculated using the MDD as the denominator and HCs as the numerator (Figure 2F).The colored bars represent the value of the density ratios.For parameter search, a = 0, the regularization parameter h varied from 0.10 to 0.01, and the Gaussian kernel width s took values of 1.2, 1.0, and 0.8.As a result, h = 0:01 and s = 1:0 were selected.After performing the density ratio estimation, it was necessary to determine the significant areas.Therefore, a permutation test was performed with a = 0, h = 0:01, and s = 1:0.The labels of HCs and MDD across all DMs were shuffled, and density ratio estimation was applied to calculate the maximum peak value 100 times (Supplementary Figure S1).Subsequently, areas above the 95th percentile of the maximum peak value were calculated (Supplementary Figures S2A, B) and clustered using DBSCAN (Supplementary Figures S2C, D).
Glass brain plots depicting the amplitude and phase of the mean DMs, histograms of frequency, and discrete-time eigenvalues within clusters in the MDD/HCs (Figure 3) and HCs/MDD (Figure 4) cases are presented.Because DMs appear in pairs with modes of identical amplitude and an anti-phase relationship, DMs at symmetric locations are paired (Figure 2E 1-2, Figure 2F 4-5, and 6-7).
In the MDD/HCs case, the glass brain plots of DM1 and DM2 were similar to those of DMN.The discrete-time eigenvalues were distributed along the unit circle, indicating stability in DM1 and DM2.The glass brain plots of DM3 were similar to those of the SMN.The discrete-time eigenvalues were relatively numerous inside the unit circle, indicating not only stability but also convergence in DM3.Additionally, because both the DMN and SMN were concurrently active in DM2, the frequency histogram was likely to show an intermediate distribution between the distributions in DM1 and DM3.
In the HCs/MDD case, the glass brain plots of DM4 and DM5 were similar to those of the VN.The discrete-time eigenvalues were distributed along the unit circle, indicating stability in DM4 and DM5.The histogram of the frequency showed a peak at approximately 0.03 Hz.The glass-brain plots of DM6 and DM7 were similar to those of the SN.The discrete-time eigenvalues were distributed along the unit circle, indicating stability in DM6 and DM7.The histogram of the frequency showed a peak at approximately 0.15 Hz.The small number of DMs in DM7 likely resulted in a negative bias of the phase and scattering of the frequency histogram.

Discussion
We devised a methodology for estimating brain-state stability across subjects by applying DMD to BOLD signals; t-SNE was applied to the DMs to disentangle the intricate curved surface spanned by the DMs into a two-dimensional manifold (Figure 2).Density ratio estimation was then performed on the twodimensional manifolds of HCs and MDD (Figures 2E, F).Consequently, it was revealed that MDD did not cause the emergence of unknown DMs from HCs but sank into specific DMs, such as DM1, DM2, and DM3.In machine learning using DMD, there are two important aspects of comparing HCs and MDD.One is interpretability in terms of physiology and the other is classification performance for biomarker.Therefore, individual-level classification between HCs and MDD was performed to demonstrate usability to the biomarker development (Supplementary Figure S6).As a result, when evaluated using 10-fold cross-validation (Supplementary Figure S7), the balanced accuracy (Bacc) was slightly better than that in the previous study (12) using sFC (Supplementary Figure S8).

Dynamic modes and cortical abnormalities of MDD
The spatial patterns of reduced DMs corresponded to the patterns observed in the cortical thickness and surface area abnormalities (29).Specifically, DM6 and DM7 exhibited spatial  Mean DM's amplitude, phase, histogram of frequency, discrete-time eigenvalue l in each HCs/MDD cluster.The left numbers correspond to the peak numbers in Figure 2. DM4 resembles a visual network (VN), has a low frequency, and is stable.DM 5 resembles a VN, has a low frequency, and is stable.DM6 resembles a salience network (SN), has a high frequency, and is stable.DM 7 resembles an SN, has a high frequency, and is stable.
map, and result in macroscopic abnormalities, such as BOLD signals (35)(36)(37)(38)(39). Related to mesoscopic phenomena, some abnormalities observed in the reuptake of neurotransmitters, such as serotonin, dopamine, norepinephrine, and GABA (40)(41)(42) resulting in neurotransmitter concentrations in plasma metabolism (43).Related to macroscopic phenomena, MDD exhibits reduced cortical thickness and surface area compared with HCs (29).As if to connect these two different scale phenomena, both the cortical abnormalities and receptor maps share similar spatial patterns (44,45).These combined abnormalities likely resulted in sinking into specific DMs, such as DM1, DM2, and DM3.Hence, if a subject transitions from HCs to MDD, it is plausible that MDD would submerge into these particular DMs alongside reductions in cortical thickness and surface area, as well as neurotransmitter reuptake abnormalities.
As a first step in integrating multiple pieces of information that reflect different aspects of MDD, it is vital to investigate the relationship between alterations in stability based on DMs and reductions in cortical thickness and surface area using large datasets.In a comprehensive study on white matter alterations in HCs and MDD, fractional anisotropy was found to be decreased in adult MDD but not significantly different in adolescent MDD compared to HCs (46).Conversely, adolescent MDD exhibited decreased cortical surface areas, particularly in regions such as the orbitofrontal cortex and lateral occipital cortex, when compared to HCs (29).Therefore, in addition to examining the structural connectivity based on the fiber structure in the white matter, it is essential to consider stability measures based on reduced cortical surface areas in both HCs and MDD.Notably, sFC can be well explained (approximately 0.9) by geometric modes (GMs) derived from the cortical geometric structure in HCs (47), suggesting that GMs could serve as a valuable stability indicator based on brain structure.
The integration of multiple indicators will be effective in psychiatric care.A combination of temporally stable trait biomarkers and temporally variable state biomarkers is necessary for early diagnosis and intervention using mechanism-based treatments (48).Therefore, structural connectivity and GMs, as temporally stable trait biomarkers, are employed as criteria for assessing stability.Additionally, DMs serve as temporally variable state biomarkers for evaluating the current cortical stability.The integration of the stability associated with cortical structural and geometric alterations and BOLD signals may shed light on previously unknown mechanisms underlying MDD.

Inconsistency with the previous studies
In MDD, negative emotions are associated with increased activity in the DMN (49) and motor impairment is associated with slow gait and slumped posture (50).Consequently, DM1 and DM2, resembling the DMN, probably emerged for experiencing negative emotions, and DM3, resembling the SMN, probably emerged for experiencing movement difficulties.
In the EL-based method (20), non-melancholic MDD tended to sink into the left CEN, whereas melancholic MDD tended to sink into both the left CEN and dorsal DMN states.In contrast, in the DMD-based method, the MDD sinks into brain states resembling the DMN and SMN.These differences can be attributed to the following three factors.First, the binarization process affects the results.In the DMD-based method, the strong amplitudes of all DMs, except for DM3, were approximately 0.03 in regions associated with the DMN, VN, and SN, and medium amplitudes were approximately 0.01 in regions associated with the SMN.In contrast, the strongest amplitudes of DM3 were associated with the SMN, but the amplitude value was only 0.003, which is approximately 1/10 compared with the other DMs.Conversely, the EL-based method requires the binarization of BOLD signals after functional network assignment to a specific region.This binarization process may have led to an outcome in which regions with amplitudes smaller than the average were considered inactive.Second, the larger number of subjects in our study may lead to more robust results than the previous study.This study included 845 subjects, whereas the previous study included 262 subjects.Lastly, regarding the subtype of MDD, this study did not differentiate between non-melancholic and melancholic MDD, whereas previous studies analyzed these subtypes separately.These methodological discrepancies and different numbers of subjects may account for the sinking into different states between the DMD-and EL-based methods.
In a large dataset study using the sFC (13), hypoconnectivities were observed within the SMN and SN, as well as between the SMN, SN, dorsal attention network (DAN), and VN in MDD.However, no significant differences were found between the DMN and frontoparietal networks (FPN).In contrast, this study identified abnormalities in the DMN, SMN, VN, and SN but no abnormalities in the DAN.A previous study using the same dataset showed that there were only a few abnormal FCs related to the DAN and many abnormal FCs related to the DMN (12).It is worth noting that the DMN and DAN exhibit an inverse correlation, wherein DAN activation leads to DMN suppression (51).Therefore, it is possible that the subjects in this study activated the DMN, while those in the larger dataset study used an sFC-activated DAN (13).

Relationships among DMs' spatial pattern, histogram of frequency, and discrete-time eigenvalue
The amplitude of DM3 exhibited a spatial pattern resembling that of the SMN and was approximately 0.003, which was approximately 1/10 smaller than the amplitudes of the other DMs.The amplitudes in DM6 and DM7 were stronger in the SN and slightly stronger in the SMN than in the other DMs.The amplitudes of DM1 and DM2 were stronger in the DMN and slightly stronger in the SMN.Consequently, the SMN tended to appear more frequently in conjunction with other networks.Furthermore, the observation that the SMN tended to co-occur with low-frequency DM1 and DM2, as well as high-frequency DM6 and DM7, suggests that DM3 transmitted information across a broad range of frequencies, resulting in a smoother frequency distribution compared to the other DMs.

FIGURE 1
FIGURE 1 FIGURE 2 Two-dimensional manifolds of HCs (A) and MDD (B) with t-SNE, kernel density estimation of HCs (C) and MDD (D), and density ratio distribution estimated by relative unconstrained least-squares importance fitting (RuLSIF) in the case of MDD/HCs (E) and HCs/MDD (F).The points on the twodimensional manifold indicate DMs (A, B).The curved lines on the density estimation indicate contour lines (C, D).The red numbers indicate the peak number.In the MDD/HCs case, the peaks located at the far left and far right were not assigned numbers due to their lack of significance at the 95% confidence level (E, F).MDD/HCs shows increased DMs in MDD, and HCs/MDD shows reduced DMs in MDD.

FIGURE 3 Mean
FIGURE 3Mean DMs' amplitude, phase, histogram of frequency, discrete-time eigenvalue l in each MDD/HCs cluster.The left numbers correspond to the peak numbers in Figure2.DM1 resembles the default mode network (DMN), has a low frequency, and is stable.DM2 resembles DMN, has a flat frequency, and is stable.DM3 resembles a sensory-motor network (SMN), has high frequency, and tends to converge.