Spontaneous transient states of fronto-temporal and default-mode networks altered by suicide attempt in major depressive disorder.

Major depressive disorder (MDD) is associated with increased suicidality, and it's still challenging to identify suicide in clinical practice. Although suicide attempt (SA) is the most relevant precursor with multiple functional abnormalities reported from neuroimaging studies, little is known about how the spontaneous transient activated patterns organize and coordinate brain networks underlying SA. Thus, we obtained resting-state magnetoencephalography data for two MDD subgroups of 44 non-suicide patients and 34 suicide-attempted patients, together with 49 matched health-controls. For the source-space signals, Hidden Markov Model (HMM) helped to capture the sub-second dynamic activity via a hidden sequence of finite number of states. Temporal parameters and spectral activation were acquired for each state and then compared between groups. Here, HMM states characterized the spatiotemporal signatures of eight networks. The activity of suicide attempters switches more frequently into the fronto-temporal network, as the time spent occupancy of fronto-temporal state is increased and interval time is decreased compared with the non-suicide patients. Moreover, these changes are significantly correlated with Nurses' Global Assessment of Suicide Risk scores. Suicide attempters also exhibit increased state-wise activations in the theta band (4-8 Hz) in the posterior default mode network centered on posterior cingulate cortex, which can't be detected in the static spectral analysis. These alternations may disturb the time allocations of cognitive control regulations and cause inflexible decision making to SA. As the better sensitivity of dynamic study in reflecting SA diathesis than the static is validated, dynamic stability could serve as a potential neuronal marker for SA.

Introduction

Depression is the leading cause of mental health burden worldwide, and in extreme cases leads to suicide [1]. A meta-analysis provides that suicide attempt (SA) has 31% lifetime prevalence in individuals with major depressive disorder (MDD) across the world [2]. The topic of depression and suicide is unprecedentedly urgent and essential as the increased suicidality in the COVID-19 pandemic reported from different countries [3-5]. At present, clinical assessments of suicide risk usually rely on retrospective information, which are often unsatisfactory [6]. Additionally, nearly 80% of patients who committed suicide didn’t report suicidal ideation in their last communication about this [7]. Therefore, there is a need to explore the potential mechanism of suicidal behavior to find suicidal biomarkers which may assist in improving the future clinical evaluation of suicidal risk [8, 9]. Identifying brain alterations by neuroimaging may help to understand the neuronal mechanism of suicidality and develop such markers [10].

Resting brain activity has been linked to higher order cognitive processes and previous studies suggested that an impairment of cognitive control may underpin the high suicide rates in individuals with MDD [11]. To better understand the neural mechanisms underlying resting-state brain functions, there have been a shift in brain mapping from the study of discrete functional regions to the identification of networks [12]. The existed evidence for resting-state network (RSN) changes in SA comes largely from functional magnetic resonance imaging (fMRI) [13-15]. For example, suicide attempters with MDD have been demonstrated with synchronous alterations in intrinsic activity of the frontal, temporal, and parietal areas [16], altered intra- and inter-network connectivity among default mode network and salience network, as well as the right frontal-parietal network [17]. Furthermore, the decreased functional network connectivity between insular–default mode network and insular–cerebellum was associated with the suicide and stress level in the suicidal depressed group [18].

Despite the above research progresses, most analysis undertaken in a static definition of the whole-time scale usually ignore the dynamic activity in resting state and erase the tiny features in particular time point. Actually, resting state has been proved to be underpinned by much richer spatiotemporal dynamics then previous studies assumed and the network could be better characterized using some extra time-varying measures of interactions [19, 20]. But to date, dynamic studies on suicide issue of MDD are still limited and they all focused on suicide ideation (SI). Dynamic whole-brain connectomics using a sliding window were found with increases in the overall topological properties among MDD patients with SI, and the features may correlate with the severity of SI [21]. Another dynamic characterization of low-frequency fluctuation (dALFF) on resting-state fMRI using sliding-window analysis suggested that the SI MDD group showed decreased brain dynamics in various regions involved in executive and emotional processing [22]. Dynamic evidence that SI in patients with MDD may be related to an abnormality in habenula was also found by sliding-window fMRI analysis [23]. These dynamic studies reported their findings on the differential regions of SI, but the instability of neuronal networks has not been considered.

In addition, the temporal resolution of sliding window approaches is limited, as each window requires relatively large amounts of data, typically several seconds in length [24]. Due to its limited temporal resolution, fMRI currently cannot be used to image dynamic brain activity in the time frame in which neuronal processes occur, i.e., in the sub-second range [25]. But it is likely that the brain supports complex thought and behavior by dynamic recruitment of whole brain networks across millisecond time-scales [24]. The signatures of these dynamics may be observable in electrophysiological recordings, which reveal that there are sub-second spontaneous states in resting state brain activity [19, 26–28]. These spontaneous transient states may organize and coordinate neural activity in brain networks. Although SA is the most important risk factor for subsequent suicide behavior, how sub-second transient states of large-scale brain networks are altered in SA is still unknown. Therefore, we require electrophysiological recordings to have more complete understanding measures of the suicidal trajectory of networks.

In the current study, we used resting-state magnetoencephalography (MEG) to track the abnormality of brain dynamics underlying SA of MDD patients. The major advantage of MEG is the ability to capture neuronal electrophysiological signals with high temporal resolution and rich frequency information. The development of analysis methods that capture dynamics, Hidden Markov Model (HMM), opens such an opportunity to explore how dynamics in oscillations across a spectral range supports the electrophysiological networks [19, 28]. HMM was developed without a prior on dynamic windows to assess network dynamics across the brain parcellations [24, 29]. To our knowledge, this HMM application is the first dynamic electrophysiological network study on SA of MDD, which could characterize the brain alternations of SA in a spatially, temporally and spectrally defined way. We hypothesized that suicide attempt in MDD may be associated with aberrant dynamics involving activation of the frontal, temporal, parietal and default mode networks. Our aim is to provide greater mechanistic insight into the brain dysfunctions by exploring the altered brain network dynamics underlying SA, and also to yield more potential measures to serve as diagnostic neuronal markers as a supplement to the clinical suicidal assessments.

Materials and methods

Participants

80 MDD patients were recruited at the Inpatient Department of Psychiatry of the Affiliated Brain Hospital of Nanjing Medical University, while 50 healthy controls (HCs) were enrolled through advertisements in the same region. Resting-state MEG scanning was performed on patients at the time of hospitalization. All subjects were right-handed and fulfilled the conditions to undergo a MEG scan. Among them, 2 MDD and 1 HC participants were excluded due to excessive head movements and poor image quality. In the final analyses, 78 MDD and 49 HC were included. Patients were assessed with the MINI-International Neuropsychiatric Interview (M.I.N.I.) [30] in accordance with the criteria of both the Diagnostic and Statistical Manual of Mental Disorders, 4th edition (DSM-IV). The severity of depression was evaluated by the 17-item Hamilton Depression Scale (HAMD-17).

MDD patients consisted of suicide-attempted (SA) subgroup including 34 participants and non-suicide (NS) subgroup including 44 participants. Patients were included in the SA subgroup when they had at least one documented self-injurious act with the intent to die in the current episode, and confirmed by medical records that the HAMD-17 third item (suicide) score ≧ 2 [6, 31]. Patients were included in the NS subgroup when they were without any suicide attempt history during the present or in previous episodes. Patients also met with the following inclusion criteria: (1) having a depressive episode with a HAMD-17 total score > 17; (2) no comorbidity with other DSM-IV axis-1 disorders such as schizophrenia, substance abuse, obsessive compulsive disorders, and generalized anxiety; (3) no use of psychotropic medication including antidepressants, antipsychotics and benzodiazepines for the past 2 weeks; (4) no physical therapy such as repetitive transcranial magnetic stimulation (rTMS) or electroconvulsive therapy (ECT) during the past six months. Participants with the following criteria were excluded from the study: (1) serious medical conditions such as organic brain disorders and severe somatic diseases; (2) history of alcohol and drug abuse; (3) pregnancy. HCs were screened for family history of any mental disorders and adopted similar exclusion criteria to that of patients.

Table 1 summarized the demographic and clinical characteristics of all participants. There is no group difference in gender, age and education years. The difference between the NS and SA subgroups in HAMD-17 total scores could be eliminated by subtracting the third (suicide) term, which implies there is no difference in their other depressive severity. In addition, the suicide risks in MDD patients were also assessed by the Nurses’ Global Assessment of Suicide Risk (NGASR) [32], which also reflects difference in suicide level between the NS and SA subgroups.

Table 1

Demography for all subjects

	MDD patients		Healthy controls	p values
	Non-suicide (NS)	Suicide-attempt (SA)
Gender (F/M)	23 M/21F	14 M/20F	25 M/24F	0.577
Age (years)	30.8 ± 8.6	28.1 ± 9.6	30.9 ± 7.2	0.287
Education (years)	13.8 ± 2.8	13.7 ± 3.0	14.7 ± 1.4	0.160
Course of disease (months)	68.8 ± 66.8	73.3 ± 70.6	−	0.429
Number of episodes of depression	3.0 ± 1.8	3.5 ± 2.8	−	0.224
Family history of mental disorder (Y/N)	13/31	10/24	−	0.990
Family history of suicide (Y/N)	2/42	2/32	−	0.791
HAMD-17 total scores	20.5 ± 4.6	23.1 ± 4.2	−	0.015*
HAMD-17 3rd item (suicide)	0.3 ± 0.4	3.1 ± 0.6	−	0.000*
NGASR scores b	6.1 ± 2.1	12.3 ± 1.9	−	0.000*

*Significant differences between groups (p < 0.05)

a Note here, the third item of HAMD-17 is to evaluate the suicide level from 0 to 4

b NGASR, Nurses’ Global Assessment of Suicide Risk

MEG image acquisition

MEG data were recorded in a magnetically shielded room using an Omega 2000, 275 channel CTF MEG system (VSM Med Tech Inc, Port Coquitlam, Canada). Resting state MEG scanning lasted for 4 min with 300 Hz sampling rate. The subjects were instructed to lie in the supine position with closed eyes and to neither fall asleep nor think of anything in particular during the scanning session.

Individual anatomical images were acquired with a Siemens Verio 3 T MRI system using a high-resolution, T1-weighted, 3D gradient-echo pulse sequence (TR = 1900 ms, TE = 2.48 ms, FA = 9°, slices number = 176, slice thickness = 1 mm, voxel size = 1 × 1 × 1 mm3, FOV = 250 × 250 mm2). Three fiducial markers, place at the nasion, left and right preauricular, enabled offline MRI and MEG co-registration.

Preprocessing and source reconstruction

The raw MEG data were preprocessed by the Fieldtrip toolbox (fieldtrip.fcdonders.nl). All-sensor signals were firstly band-stop filtered to remove power-line inference (50 Hz), and then band-pass filtered to 1–75 Hz frequency range. Within fieldtrip interface, trials and channels deviating in the overall distributions were removed. Further step was applied using temporal Independent Components Analysis (tICA) across the sensors using the FastICA algorithm. Artifacts related to breath, heart beat and muscle movement were rejected by visual check in this step. Then a linearly constrained minimum variance (LCMV) beamformer was utilized to project the resulting pre-processed MEG data onto a regular 8 mm grid source space [33].

Parcellation and orthogonalization

The brain was parcellated with a weighted parcellation of 39 cortical regions, which was identified from a resting-state ICA decomposition of fMRI data from the first 200 subjects of the Human Connectome Project and had been previously used in various MEG studies [24, 34, 35], particularly in 275-channel CTF MEG data [36]. Single time series as the first principal component was acquired to represent each parcellation. A multivariate symmetric orthogonalization was then adopted to attenuate the spatial leakage effects [36].

The hidden Markov model

The HMM assumes that a time series can be described using a hidden sequence of a finite number of states [37]. At each time point, only one state is active, the probability of a state being active at time point t is modelled to be dependent on which state was active at time point t − 1 (i.e. it is order-one Markovian) [29]. HMM aims to discover these hidden brain states as well as the likely sequence of transitions between them. The HMM could be inferred from the source-space MEG data using the HMM-MAR toolbox (https://github.com/OHBA-analysis/HMM-MAR). To alleviate overfitting issues, time delay embedded HMM (TDE-HMM) [28, 29] applies a different variety of the HMM to time courses and infers a multivariate Gaussian distribution describing a delay-embedding of the source time courses.

HMM inference

The source-reconstructed time courses were embedded with time lags varying between -7 and 7. In our dataset with 300 Hz sampling rate, it corresponded to windows of 100 ms and specified a 50 ms lag in both directions. The lags resulted in (N parcels N lags) N timepoints data matrix for each subject. N is referred to 39 parcels here. Afterwards, the first dimension of this matrix containing the spatial and lag information was reduced by projecting the matrix onto the first 4 N parcels components of a principal component analysis. The HMM-MAR uses stochastic inference, which is based on taking subsets or batches of subjects at each iteration instead of the entire data set [38]. The HMM-MAR toolbox uses the alternative to the standard HMM, a stochastic variational inference approach, that can be applied to neuroimaging datasets, by greatly reducing its computational cost [39]. The batch size was set to use 15 continuous data segments at each iteration and number of variational inference cycles was set to 500. The whole HMM framework ran 10 times to ensure the stability of results and the best performance with lowest free energy was accepted here.

HMM temporal parameters

After HMM analysis, time series of posterior probabilities was inferred to represent the occurrence probability of a state at a time point. After the inference process, the Viterbi path was also computed [40]. This is defined as the most probable sequence of states representing the mutually exclusive state allocations. They were used to calculate the temporal parameters of each state, includes Fractional Occupancy (FO), Life Time (LT) and Interval Times (IT). These parameters could provide an overview of brain state dynamics captured by HMM [19]. FO refers to the proportion of each state time spent in the total time length. LT is the mean dwell time of each state on a single visit and IT is the mean time between the state visits [12]. In the current step, the above temporal parameters of each HMM state were calculated for each subject individually.

HMM state-wise spectral analysis

To acquire power for each frequency band and HMM state, we made use of a nonparametric estimation, using a novel state-wise version of the multi-taper [28]. The state-wise power across the parcellations were computed from the results of HMM in the spectral range of broadband (1–30 Hz) for each state. The spectral content and spatial maps were acquired for each subject and then averaged among all subjects to see the correspondence of HMM states to intrinsic brain networks. Afterwards, we factorized the spectral information in each subject into different frequency modes for ease of group comparison in each frequency band.

Statistical analysis

The overall dynamics of three temporal parameters (FO, LT and IT) were preliminarily compared between the whole MDD cohort and the HC cohort, and then between the SA subgroup and NS subgroup of MDD cohort to find suicidal differences. These differences [HC vs. MDD, MDD(SA) vs. MDD(NS)] were all found with permutation tests and the resulted p values were false discovery rate (FDR) corrected. Pearson correlation analysis was then used to study the associations among suicide risk measures (NGASR scores) and temporal parameters of different states in SA individuals.

Similar to the dynamic parameters, the state-wise power in various frequency bands (theta, alpha, beta) were also compared between cohorts [HC vs. MDD, MDD(SA) vs. MDD(NS)] with permutation test and followed by FDR correction.

Static spectral analysis

As a comparison to the above dynamic methods, the static frequency analysis was calculated with the same multi-taper method. Static power spectra were calculated directly from the whole-length time courses of all 39 regions without any state information. The static power spectra were compared between cohorts [HC vs. MDD, MDD(SA) vs. MDD(NS)] in the exact same way with the dynamic state-wise spectra.

Results

Eight HMM states were identified by TDE-HMM inference on all subjects. The mean activation maps and temporal parameters of eight HMM states across all subjects are displayed in Supplementary Figure S1. The goal of the variational inference in the current HMM analysis is the minimization of the so-called free energy. Free energy from the model inference provides an approximation to model evidence (accuracy—complexity) [39]. The percentage change decreases as the number of states increases and the improvement becomes negligible beyond eight states as shown in the supplementary Figure S2. In addition to the objective measure, we also followed the approach in [19] to determine the final HMM states (K) by testing a range of values for K. The corresponding results for mean activation maps are shown in supplementary Figure S3. The principle here is not to establish the ‘correct’ number of states, but to identify an optimal number that provides a reasonable description of the dataset for each specific purpose [24]. Our selection of state number as eight here is also driven by that the spatial topographies are in consistent with previous MEG HMM studies [12, 19] and traditional resting state networks.

Changed temporal dynamics in suicide attempters

We found differences in temporal parameters of state 1, 4, 7 and 8 between the whole MDD cohort and HC cohort (Supplementary Figure S4). Then we further explored differences between SA and NS subgroups of MDD to find the suicidal features. The changed temporal dynamics related with SA were found in state 4 (Fronto-temporal), 7 (Sensorimotor) and 8 (Parietal) as present in Fig. 1. For the fronto-temporal network, the FO of SA patients is significantly higher than that of the NS patients (p < 0.001 after correction), and IT is significantly lower (p < 0.001 after correction). The significant higher FO and LT also are displayed for the parietal network in SA patients accompanying with lower IT (p < 0.001 after correction). And for the sensorimotor network, LT is lower in the SA patients (p < 0.001 after correction).

Fig. 1

A The spatial maps for state 4 (Fronto-temporal), 7 (Sensorimotor) and 8 (Parietal) for all subjects. Please note that the activation maps have been thresholded here to visualize. B Corresponding plots for each state show fractional occupancy, state lifetimes and interval times between suicide attempted (SA) and non-suicide (NS) MDD subgroups. Asterisks (***) denote significantly changed temporal dynamics with p < 0.001. The crosses in the figure represent the mean values

Furthermore, the temporal parameters of fronto-temporal network (state 4) in SA subgroup of MDD could be correlated with their suicide risks, which were clinically measured by NGASR. As shown in Fig. 2, FO of the fronto-temporal network is positively correlated with suicidal risk (r = 0.4239, p = 0.0125) and IT is negatively correlated (r = -0.4384, p = 0.0095). The FO and IT are validated to be not correlated with severity of depression within SA patients.

Fig. 2

Significant correlations between suicide risk scores (NGASR) and dynamic parameters of fronto-temporal network (fractional occupancy and interval times) in suicide attempted MDD patients. NGASR, Nurses’ Global Assessment of Suicide Risk

Changed state-wise power in suicide attempters

Significant differences on state-wise power were found in state 1, default mode network (DMN) between the whole MDD cohort and HC cohort as shown in Supplementary Figure S5. To detect suicidal features, spectral activation differences of the DMN between SA and NS subgroups of MDD were further compared in each independent frequency band. DMN (state1) shows an abnormal power activation pattern in theta band (Fig. 3A) in SA patients compared with NS patients. In the posterior part of DMN, SA patients activate with higher power in theta band than NS depressed patients do (after correction). Interestingly, the posterior cingulate cortex (PCC), on which the posterior DMN centers [41], is included in the differential regional sets (p < 0.001 after correction). In Fig. 3B, power of the whole brain averaged among 39 parcellations are plotted for 8 dynamic states (acquired from HMM inference) and a static condition (calculated in the whole-length time scale without HMM). The power of PCC is also plotted in the same way. The power of state 1 (DMN) is mainly activated in the low frequency and much higher than all the other states here. The same situation is applicable to the PCC. The PCC power is also highest in state 1 compared with all the other states and the static condition. Suicide induced spectral changes in the static and HMM-specific frequency content are displayed in Fig. 3C. The dynamic spectral analysis is likely to decompose the static spectral findings in different temporal HMM states, but could reveal more clear and detailed differences. To validate that, we also present the static spectral analysis in the following part.

Fig. 3

A Comparison between suicide attempted and non-suicide MDD subgroups on the power of default mode network derived from HMM states. Colorbars represent -log10 transform of p values, which mean 1.3–3.5 in the colorbars correspond to p values in 0.05–0.0003. B Powers of the whole brain averaged across regions for the 8 states in HMM dynamic analysis and for the whole time-scale static analysis. Power distributions for the PCC are also plotted in the same way. C Power changes in the HMM-specific and static frequency content induced by suicide attempt. The significant difference of default mode network (state1) in Figure A is boxed in red

Static spectral findings

As a comparison, static spectral power was calculated and compared also in three frequency bands, which could be found in Fig. 4A. Among the 39 parcels, the main difference between the whole MDD group and HC group is located in the prefrontal power (after correction) of theta and beta. Then the static powers of three frequency bands were further compared between MDD subgroups of SA and NS as shown in Fig. 4B. Although certain changing trends could be found in the frontal and sensorimotor areas between groups, there is no significant result survived after multi-comparison correction. The static method is not as sensitive as the HMM method to find the distinguishing characteristics of suicide attempt.

Fig. 4

A Static spectral comparisons between the whole MDD and HC cohorts. Significant decreased theta power and increased beta power in the depressed group are displayed. Colorbars represent -log10 transform of p values, which mean 1.3–3.5 in the colorbars correspond to p values in 0.05–0.0003. B Static spectral comparisons between the non-suicide and suicide attempted MDD subgroups. Changing trends could be found between MDD patients without suicide and with suicide attempt, but no significant finding could be found after correction. Colorbars represent -log10 transform of p values, which mean 1.3–2 in the colorbars correspond to p values in 0.05–0.01

Discussion

The study investigated and compared dynamic functional networks in MDD patients with and without suicide attempts. By the comparison between SA and NS subgroups of MDD, we show that FO of fronto-temporal network is a dynamic variable increasing significantly in SA patient subgroup. Furthermore, we show that this variable is a sensitive marker for SA, as the fronto-temporal FO is positively correlated with the clinically assessed suicide risk. Also in the SA subgroup, the DMN exhibits a higher state-wise activation in theta band compared with the NS group, which could not be captured by the static analysis and emphasizes the essentiality of dynamic analysis. The current study provide a more comprehensive insight into neural activities, which may be ignored via averaging those activities over time in static analysis.

A recent review has also proposed that the suicidal brain across psychiatric diagnoses seems to heavily involve the dysfunction of the fronto-temporal network [42]. Such dysfunction has also been further enriched from the dynamic aspect from our findings. In this study, dynamic parameter of the fronto-temporal network, particularly the FO, is significantly increased in the SA subgroup and positively correlated with the increases in suicide risks. From the aspect of cognitive impairment, suicide attempters were proven with impaired executive function, attention, and memory, which was linked to prefrontal lobe dysfunction [18]. The specific brain network of fronto-temporal regions was suggested for processing emotional prosody [43]. In addition, planning and acting out suicidal impulses in response to mental pain were also associated with increased activity in the frontal cortex [44]. It’s interesting that we report the increase of fronto-temporal FO in SA patients is accompanied by a decrease in IT, but with no apparent difference in state LT between the two MDD subgroups. Once the fronto-temporal state occurs, its dwell time is not changed with SA. However, the increased FO and decreased IT imply that the resting-state activity of SA patients switches more frequently into the fronto-temporal state, which may disturb the regular time allocation of the whole cognitive processing and increase mental pain in MDD patients. The correlation between temporal parameters and suicidal risk also implies further that FO of fronto-temporal network could be a potential biomaker to predict the risk level of suicide behavior.

In addition to the frontal and temporal cortices, it is agreed that dysfunction associated with the parietal cortices is also implicated in the suicidal brain. Fronto-parietal network is involved as a part of cognitive control network which is critical for problem-solving and executive functioning [14]. Previous in vivo and post-mortem neurobiological studies have also implicated ventromedial and parietal regions in suicide risk [45]. According to our HMM analysis, three temporal parameters of the parietal network are all significantly changed by SA. Compared with NS MDD subgroup, SA MDD subgroup exhibits increased FO and LT. There may be major dynamic disturbance to the related cognitive control function, so as to increase the vulnerability to suicide for MDD patients. In addition, this study found decreased LT of sensorimotor network in the SA subgroup. Sensorimotor-related areas involve the lateral temporal and occipital lobes, and primary sensorimotor cortex [46]. Attention and other cognitive processes also play a role in the modulation of rhythmic activity in sensorimotor regions and similarly may reflect sensory gating mechanisms involved in motor preparation or anticipation of sensory input [47].Based on the suicidal changing trend in the parietal and sensorimotor state, we propose that sensorimotor network works as a complement to the parietal network in the suicidal cognition processing, and thus, it exhibits the opposite changing dynamics to the parietal abnormalities in SA patients.

Of our HMM analysis, state 1 exhibits significant spatial similarity with the stationary activation pattern. Its spatial morphology, coupled with low FO and high IT (Supplementary Figure S1) between state visits suggests that this state is equivalent to what Baker and colleagues [19] termed the DMN. In addition, its power spectra is also much higher than all the other states as shown in Fig. 3B, which is in agreement with Diego et al.[29]. Our dynamic analysis finds SA related theta-band power difference in the posterior activations of DMN (Fig. 3A), including its centered region PCC. The DMN plays an important role in “internal mentation”—the introspective and adaptive mental activities in which humans spontaneously and deliberately engage [48]. Several lines of evidences suggest a perturbed sense of self in people who attempt suicide [49, 50]. Suicidal activity in mood disorder patients may be also a consequence of impaired self-referential thought processing [15]. In particular, posterior DMN has been implicated in consciousness and memory processing [41], which are both critical to suicidal susceptibility. Our results in the posterior part of DMN, SA patients activate with higher power in theta band. The evidence from EEG study reported that theta power in the fronto-central region was significantly increased in the high SI group [51]. Although the EEG study were conducted in the sensor level and reported scalp electrode changes, it supports the increased activity in theta band has clinical potential as a biomarker for identifying suicide. Thus, we validate the essentiality of exaggerated posterior DMN activation in theta band underlying SA.

To our knowledge, only few previously published studies examined suicide in MDD using MEG. Optimal use of the rich information content of MEG signals often benefits from reconstruction of the generators of the signals [52]. The previous analysis from our group of resting-state MEG data in source space found caudothalamic coupling abnormality in the high risk suicidal group [53]. In a whole-brain source-level MEG study, the gamma power of anterior insula regions was found to be associated with SI and Ketamine could improve the clinical symptoms of SI [54]. Healthy controls were not included. In the recent MEG study, they didn’t capture the group difference, but revealed that electrophysiological connectivity, as coupling of low frequency power (delta and theta) with alpha and beta power, was strongly related to SI [55]. These published MEG study on suicide didn’t consider the sub-second dynamics in spontaneous fluctuations, but validated the potential of MEG to explore SA and inspired us to make full use of MEG temporal information. As we adopted the difference analysis methods from the previous MEG application, the results are hard to be comparable with them. In our static analysis shown in Fig. 4, the difference of brain regional power can only be found between the whole MDD and HC cohorts. No static power difference between SA and NS subgroups is significant after correction. That is not surprising because the difference between subgroups in a cohort is usually not as significant as that between two cohorts. Dynamic differences in the time scale may be erased in the traditional static analysis. This dynamic analysis complements the past static findings and provides more sensitive measures. The HMM is a probabilistic model, which means it does not specify a model fit with Gaussian noise distributed around it, instead the state-specific multivariate Gaussian distribution is itself the actual model description. So we can look at the free energy, as an approximation of the model evidence, i.e. prob (model | data) [39]. As shown in supplementary Figure S2, the improvement in free energy over the K = 1 model as evidence that it is doing a better job than a static model at least. Together with the process of HMM validation completed in the supplementary material, we could be confident that the HMM was focused on relevant dynamics of MEG.

In the current study, we used the mean temporal parameters within each subject in the statistical analysis following the previous resting-state neuroimaging HMM studies [12, 56–59]. Before that, each individual does not only generate the single value of LT or IT from the direct HMM output other than FO. These application studies of HMM could support the feasibility of within-subject averaging parameters in our group comparisons, and the high dependency of the values from the same subject could be eliminated in this way. However, for analysis of task data, the scanning time of each subject may be divided into different period according to the stimulus onset and response time. In this case it will be tricky to average the dynamic parameters in different visits to get the within-subject mean values. To overcome the high dependency across the values from the same individual, a mixed-effect modeling/regression and controlling some predictors in the regression would be good and essential for analysis of task data where the task design may limit the summary statistics of state sequences. Controlling some predictors, like demographical variables could also be discussed in our current study. As noted in Table 1, demographical variables were not significantly different between groups in our dataset. Furthermore, we repeated the between-group comparisons of dynamic parameters after regressing out age and education years and performing inference on the resulting residuals. Between-group differences were unchanged when inference was performed on the residuals, suggesting that these matched variables did not significantly influence our findings. We would like to recommend readers to consider these points in their following work in such cases.

Some limitations need to be discussed for the current study. Firstly, our sample size is relatively small due to the difficulties in MEG acquisition. But it’s relatively considerate comparing with other MEG studies. Secondly, we only performed a cross-sectional study, and a longitudinal study is also in need both with neuroimaging recordings and clinical follow-up to ascertain the predictions of SA in the future. A bigger dataset with more clinical data would also be of interest in the future suicidal study.

Conclusion

Overall, we propose that fast transient dynamics is potential for capturing changes in state-wise temporal-spatial-spectral patterns, which may lead to dysfunctional emotional processing and cognitive control in SA of MDD individuals. Dynamic measures (FO and IT) of fronto-temporal network may become potential biomarkers of SA in MDD patients, and state-wise activation of DMN in specific frequency band should also be emphasized in the future dynamic studies on SA. The whole framework could fuel our ability to understand the neuronal mechanism of SA, and potentially response to the clinical need of suicidal prediction in the MDD population.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 1627 KB)