Investigating the predictive models of efficacy of accelerated neuronavigation-guided rTMS for suicidal depression based on multimodal large-scale brain networks

Introduction

Depression is a prevalent psychiatric disorder associated with high suicide risk and substantial illness burden. According to the World Health Organization (WHO), approximately 280 million people worldwide currently suffer from depression and >700,000 people die by suicide each year. Suicide ranks as the fourth leading cause of death among individuals aged 15–29 years (WHO). Accelerated neuronavigation-guided high-dose repetitive transcranial magnetic stimulation (NH-rTMS) has recently been shown to effectively and rapidly reduce suicidal ideation and alleviate depressive symptoms (Pan et al., 2020). However, the predictive models of treatment response to the accelerated NH-rTMS need to be further investigated, which is crucial for screening out MDD patients who are most likely to benefit from accelerated NH-rTMS treatment.

Despite many positive clinical findings, the overall efficacy of rTMS for major depressive disorder (MDD) remained uncertain. A previous study showed that only about 40 % of patients with depression responded to rTMS (Fitzgerald et al., 2016) with the response varying across clinical and demographic factors. However, none of these factors strongly justified patient exclusion from rTMS treatment. Benefiting from advances in neuroimaging technology, more and more studies have explored rTMS-related brain biomarkers to help accurately guide treatment decisions for MDD patients and identify those suitable for rTMS. Avissar et al. found that higher functional synchronization between the left dorsolateral prefrontal cortex and striatum at baseline was correlated with greater reductions in 24-item Hamilton Depression Rating Scale (HAMD-24) after 5 weeks of daily rTMS treatment (Avissar et al., 2017). Salomons et al. reported correlations between cortico-thalamic-striatal functional synchronization at baseline and reductions of depressive symptoms after 4 weeks of rTMS treatment (Salomons et al., 2014). These findings indicate the ability of baseline brain connectivity for predicting the antidepressant effect of rTMS. Despite of these advances, there are no studies so far that have explored neuroimaging biomarkers for predicting the treatment response to accelerated NH-rTMS in MDD patients.

The human brain is a complex system that can be modeled as a network of interconnected brain regions (Bullmore & Sporns, 2009; He & Evans, 2010). This network exhibits an intrinsically modular organization, where the brain can be parcellated into several coherent and distributed subnetworks that balance integration and segregation among brain regions (Betzel & Bassett, 2017). Previous studies have shown that different subnetworks exhibited unique topological organization that cannot be fully captured by that of the whole-brain network (He et al., 2009). In MDD, different subnetworks are also demonstrated to show specific susceptibilities with respect to their connectivity (Ghanbari et al., 2022; Javaheripour et al., 2021; Kaiser et al., 2015). For example, a meta-analysis of resting-state functional synchronization showed that MDD was mainly characterized by hypoconnectivity within the frontoparietal network (FPN) and between the FPN and dorsal attention network (DAN) but hyperconnectivity within the default-mode network (DMN) and between frontoparietal control systems and regions of the DMN (Kaiser et al., 2015). Based on these previous findings, in this study we hypothesize that different subnetworks would be associated with different capabilities in predicting the treatment response to accelerated NH-rTMS in MDD patients. Constructing predictive models at the subnetwork level can provide more granular insights into how baseline brain networks are associated with the treatment efficiency of accelerated NH-rTMS in MDD patients.

The aim of this study was to establish predictive models of treatment response to accelerated NH-rTMS in MDD patients with suicidal ideation based on multimodal large-scale brain networks at the subnetwork level. To this end, we first collected multimodal magnetic resonance imaging data from 32 treatment-naive MDD patients with suicidal ideation before accelerated NH-rTMS treatment. Treatment outcomes were defined as the changes of HAMD-24 and Beck Scale of Suicidal Ideation (BSSI) scores after 7-day accelerated NH-rTMS treatment. Then, individual morphological, structural, and functional brain networks were constructed, for which brain connectivity were extracted within each and between each pair of seven putative subnetworks. Finally, linear support vector regression methods were utilized to examine the ability of subnetwork-level brain connectivity in predicting antidepressant and anti-suicidal effects of accelerated NH-rTMS.

Methods

Participants

Thirty-two treatment-naive MDD patients with suicidal ideation were recruited in this study. Written informed consent was obtained from all participants. Each participant and their parents or legal guardians were provided with a comprehensive overview of the study. All the patients fulfilled the following criteria: 1) aged between 13 and 45 years; 2) the Han ethnicity; 3) right-handedness; 4) met the Diagnostic and Statistical Manual of Mental Disorder (DSM) criteria for current unipolar MDD, as confirmed by the Mini-Neuropsychiatric Interview (MINI) conducted by two qualified psychiatrists (Sheehan et al., 1998); 5) the HAMD-24 score higher than 20 (Hamilton, 1967); 6) the Beck Scale of Suicidal Ideation (BSSI) score higher than 12 (Beck et al., 1979); and 7) the score for Item #3 of the HAMD-24, which pertains to "suicide", higher than 3. The exclusion criteria for the MDD patients were as follows: 1) any treatment received prior to the study; 2) history of seizures, head trauma or other brain abnormalities; 3) significant medical illness; 4) psychotic symptoms, alcohol or substance abuse; 5) any current or past psychiatric axis-I/-II disorders (excluding MDD); 6) contraindications to MRI scan, including metallic implants, retractors or braces, and claustrophobia; and 7) women who were pregnant or lactating.

Accelerated NH-rTMS procedure

The accelerated NH-rTMS procedure targeting Brodmann Area (BA) 46 (MNI coordinates: −44, 40, 29) was performed in accordance with our previous studies (Pan et al., 2023, 2020; Zhao et al., 2024). Briefly, each individual MRI scan was used to generate a 3D curvilinear reconstruction of the brain, which was processed using Brainsight software (version 2.3; Rogue Research Inc., Montreal, QC, Canada) within the Brainsight TMS navigation system. Treatment was administered by a qualified practitioner independent of the outcome evaluator, using a Magstim Rapid² rTMS device (Magstim Company Ltd., Whitland, UK) equipped with an eight-figured coil. Specifically, the accelerated NH-rTMS regimen was performed over seven consecutive days. Each daily session consisted of 120 trains (5 s per train, 10 Hz) with 15-second inter-train intervals, totaling 6000 pulses per session. During the treatment, stimulation intensity was set at 100 % of the resting motor threshold (rMT), adjusted based on participant tolerance. The rMT was defined as the minimum stimulation intensity required to elicit a motor response in the contralateral abductor pollicis brevis muscle at rest.

Assessments

The severity of suicidal ideation and depression was assessed at baseline and day 7 post-treatment using the BSSI and HAMD-24, respectively. The percentage changes of total BSSI and HAMD-24 scores were further calculated to measure treatment outcomes.

MRI data acquisition

MRI data of the MDD patients were collected at baseline, including T1-weighted structural imaging, diffusion tensor imaging (DTI) and resting-state functional magnetic resonance imaging (rs-fMRI). All the MRI data were acquired with the Signa HDxt 3.0T MRI system (GE Healthcare, Milwaukee, WI). The high-resolution T1-weighted images were acquired using a 3D sagittal BRAVO sequence with the following imaging parameters: repetition time (TR) = 7.8 ms, echo time (TE) = 3 ms, flip angle = 7°, field of view (FOV) = 256 × 256 mm², slice thickness = 1 mm, voxel size = 1 × 1 × 1 mm³, and 192 slices. The DTI images were acquired using a single-shot SE-EPI sequence with the following parameters: TR = 13,000 ms, TE = 86.4 ms, flip angle = 7°, FOV = 256 × 256 mm², slice thickness = 3 mm, voxel size = 2 × 2 × 3 mm³, 48 slices, diffusion gradient directions = 32, b-value = 1000 s/mm², and 3 unweighted scans (i.e., b-value = 0). The rs-fMRI images were acquired using a 2D gradient echo EPI sequence with the following imaging parameters: TR = 2000 ms, TE = 30 ms, flip angle = 90°, FOV = 220 × 220 mm², slice thickness = 4 mm, voxel size = 3.44 × 3.44 × 4 mm³, and 33 slices.

MRI data preprocessing

T1-weighted images. All structural images were preprocessed using the Computation Anatomy Toolbox 12 (CAT12, version r1113) based on the Statistical Parametric Mapping 12 (SPM12, version 6685). The CAT12 provides a volume-based method to estimate cerebral surface morphology, eliminating the need for extensive cortical surface reconstruction and thus saving time. First, the structural images were segmented to obtain white matter, gray matter, and cerebrospinal fluid. The gray matter maps were then normalized to the standard Montreal Neurological Institute (MNI) space and subsequently modulated to correct for changes in tissue volume during the normalization. Finally, resampling was performed for resultant gray matter volume (GMV) images with a resolution of 1.5 × 1.5 × 1.5 mm³. In addition to GMV, we computed four surface-based morphological indices, including cortical thickness (CT), fractal dimension (FD), gyrification index (GI), sulcal depth (SD). After the tissue segment, the projection-based thickness method was applied to estimate CT (Dahnke et al., 2013), and a central surface was created. A topology correction using spherical harmonics was applied to the derived central surface to address topological defects (Yotter, Dahnke et al., 2011). By using spherical harmonic reconstruction, FD, GI, and SD were calculated as the slope of a logarithmic plot of surface area versus the maximum l-value (Yotter, Nenadic et al., 2011), the absolute mean curvature (Luders et al., 2006), and the Euclidean distance between the central surface and its convex hull, respectively. Subsequently, individual morphological maps were resampled to the fsaverage template and smoothed with a Gaussian kernel (15-mm full width at half maximum for CT and 25-mm for FD, GI, and SD). The larger smoothing kernel sizes for FD, GI, and SD reflects nature of these folding measures, which incorporate contributions from both sulci and gyri, necessitating a filter size that exceeds the distance between a gyral crown and a sulcal fundus.

DTI. The DTI data were preproceed with the Pipeline for Analyzing brain Diffusion images (PANDA) toolbox (Cui et al., 2013) based on FMRIB's Software Library (Smith et al., 2004). Specifically, the preprocessing consisted of the following steps: 1) conversion of DICOM files to NIfTI format; 2) brain mask generation; 3) cropping of the raw images by removing non-brain areas; 4) co-registration of individual diffusion-weighted images to their corresponding mean b0 image using an affine transformation in order to correct for head motion and eddy current-induced distortions; and 5) estimation of the diffusion tensor model at each voxel followed by diagonalization to obtain three eigenvalues and corresponding eigenvectors (Basser & Pierpaoli, 1996).

Functional images. The rs-fMRI data were preprocessed with the GRETNA toolbox (J. Wang et al., 2015) based on the SPM12 package. Specifically, the preprocessing consisted of the following steps: 1) deletion of the first 10 vol to allow for signal stabilization; 2) correction for slice-timing differences and head motion; 3) spatial normalization to the MNI space based on the transformation fields derived from tissue segmentation of T1-weighted images; 4) resampling to a 3-mm isotropic resolution; 5) linear detrending; 6) band-pass filtering (0.01–0.08 Hz) to reduce the effects of low-frequency drift and high-frequency physiological noise; 7) removal of the effects of nuisance variables on voxel-wise time series, including white matter signals, cerebrospinal fluid signals, global mean signals, and 24-parameter head-motion profiles. Four participants were excluded due to displacement larger than 2 mm, rotation larger than 2.5°, or a mean framewise displacement larger than 0.2 mm. Notably, the nuisance signals were filtered using the same frequency intervals as for the step 6 to prevent reintroduction of non-interest frequencies and improve noise correction within the relevant frequency range (Hallquist et al., 2013).

Networks construction

Cortical parcellation. To define brain regions, we utilized the Schaefer2018 atlas to parcellate the cortex surface into 200 regions of interest (ROIs), with each ROI as a node in the networks (Schaefer et al., 2018). These ROIs can be divided into seven subnetworks: visual network (VN), somatomotor network (SMN), DAN, ventral attention network (VAN), limbic network (LN), FPN and DMN (Yeo et al., 2011).

Morphological similarity. Five morphological brain networks were constructed for each structural image, termed as GMV-, CT-, FD-, GI-, and SD-based networks. As in our previous studies (Li et al., 2021, 2023; Wang et al., 2016), a Jensen-Shannon divergence (JSD)-based approach was then applied to measure interregional morphological similarity. Specifically, values for each morphological index within each ROI were extracted to estimate a probability density function using a MATLAB function, ksdensity. After the conversion of the estimated probability density function to probability distributions (PDs), the JSD was calculated between regional pairs based on their PDs. Formally, for two PDs P and Q, the JSD was calculated as:

where KLD (P∥Q) = ∑i=1nP(i)logP(i)Q(i) with n being the number of sampling points (2⁸ in this study). The JSD was further converted to a distance metric (square root) and subtracted it from 1, producing values that indicated similarity between two PDs, with 0 representing complete different and 1 representing exact similar.

Structural connectivity. Three structural brain networks were constructed for each participant, termed as fiber number-, fractional anisotropy (FA)-, and fiber length-based networks. Specifically, whole-brain white matter fiber tracts were reconstructed from the preprocessed DTI data using a continuous streamline-tracking algorithm (Mori et al., 1999). The tractography was terminated when it reached a voxel with a FA <0.2 or when the turning angle exceeded 45° between adjacent steps. Two ROIs were considered to be structurally connected if at least one streamline had endpoints located in both ROIs. The number, averaged FA, and averaged length of fiber tracts between each pair of ROIs were extracted to quantify the corresponding structural connectivity. Notably, the Schaefer2018 atlas was warped from the MNI space to individual native diffusion spaces where whole-brain white matter fiber tracts were reconstructed. This was done via the inverse transformation of registration from individual T1-weighted images to their b0 images (linear) and then to the MNI space (nonlinear) (Gong et al., 2009).

Functional synchronization. For each ROI in the Schaefer2018 atlas, the mean time series was first extracted. Then, the Pearson correlation was calculated to quantify the functional synchronization between each pair of regional mean time series, resulting a 200 × 200 functional brain network for each participant.

Prediction of treatment outcomes

Linear support vector regression methods were utilized to predict treatment outcomes of MDD patients (percentage changes of total BSSI and HAMD-24 scores) at the subnetwork level using the GMV-, CT-, FD-, GI-, and SD-based morphological similarity, fiber number-, FA-, and fiber length-based structural connectivity, and functional synchronization, respectively. The brain connectivity within each subnetwork and between each pair of subnetworks was used for the prediction, respectively. The predicted values of treatment outcomes were calculated via a leave-one-out cross-validation procedure and the model performance was evaluated by right-tailed Pearson's r between the real treatment outcomes and predicted values. A false discovery rate (FDR) procedure was used to correct for multiple correlations across 28 predictions (7 within-subnetwork brain connectivity and 21 between-subnetwork brain connectivity) for each type of brain networks at the level of q < 0.05. Notably, before constructing the predictive models, only brain connectivity showing significant correlations (Spearman correlation, uncorrected P < 0.05) with individual treatment outcomes in the training set were used.

To find out which regions contributed the most to the prediction, we first extracted the brain connectivity that was consistently used across all training sets during cross-validation. Then, the linear coefficients of these connectivity derived from the linear support vector regression model were averaged across all the training sets to quantify the contributions of brain connectivity to the prediction. These contributions were further transformed to Z-scores. Finally, for a given region, its contribution to the prediction was calculated as the sum of the absolute value of contributions across all connections linked to the region.

To examine the effects of imaging modality on the prediction performance (r value), one-way analysis of variance was used with morphological similarity, structural connectivity, and functional synchronization as three levels of the main factor. Furthermore, we examined whether specific subnetworks outperformed the other subnetworks for the prediction. Specifically, for a given subnetwork, the prediction performance (r value) of the models constructed by the brain connectivity within the subnetwork and between the subnetwork and each of the other subnetworks were compared with the prediction performance (r value) of the models constructed by the brain connectivity that were not related to the subnetwork using a right-tailed two-sample T test. The significance levels of all comparisons were estimated using permutation tests (10,000 times). The FDR procedure was used to correct for multiple comparisons for examining the effects of imaging modality (3 comparisons) and for examining the effects of subnetwork (7 comparisons).

Results

Prediction for changes of total BSSI scores

Morphological similarity. The GMV-based morphological similarities between the VN and DAN (r = 0.464, P = 0.004), between the VN and LN (r = 0.512, P = 0.002), and between the DAN and FPN (r = 0.518, P = 0.001), the CT-based morphological similarities between the SMN and DAN (r = 0.502, P = 0.002) and between the DAN and DMN (r = 0.610, P < 0.001), and the FD-based morphological similarities within the FPN (r = 0.575, P < 0.001) significantly predicted the percentage changes of total BSSI scores (Fig. 1). Among the GMV-, CT-, and FD-based predictions, the mostly contributed regions were the left middle occipital gyrus in the VN (Z = 4.585), the right frontal eye fields in the DAN (Z = 6.260), and the left lateral prefrontal cortex in the FPN (Z = 1.829), respectively (Fig. 2). No significant predictions were observed for the GI- or SD-based morphological similarities (P > 0.05, FDR corrected).

Fig. 1.

Prediction performance for the changes of total Beck Scale of Suicidal Ideation scores. *, significant linear correlation between real scores and predicted values (P < 0.05, FDR corrected). CT, cortical thickness; DAN, dorsal attention network; DMN, default mode network; FA, fractional anisotropy; FD, fractal dimension; FPN, frontoparietal control network; GI, gyrification index; GMV, gray matter volume; LN, limbic network; rs-fMRI, resting state-functional magnetic resonance imaging; SD, sulcal depth; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.

Fig. 2.

Connectional contributions to the predictions for the changes of total Beck Scale of Suicidal Ideation scores. Regional sizes are proportional to total absolute contributions across the linked connectivity. CT, cortical thickness; DAN, dorsal attention network; DMN, default mode network; FD, fractal dimension; FPN, frontoparietal control network; GMV, gray matter volume; LN, limbic network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.

Structural connectivity. The fiber length-based structural connectivity between the SMN and DAN (r = 0.633, P < 0.001) significantly predicted the percentage changes of total BSSI scores (Fig. 1). The left rolandic operculum and insula in the SMN (Z = 2.473) contributed the most to the prediction (Fig. 2). No significant predictions were observed for the fiber number- or FA-based structural connectivity (P > 0.05, FDR corrected).

Functional synchronization. No significant predictions were observed for the functional synchronization (P > 0.05, FDR corrected) (Fig. 1).

Effects of imaging modality. Significant effects of imaging modality were observed on the model performance of predicting the percentage changes of total BSSI scores (F = 3.9657, P = 0.022), with morphological similarity (T = 2.563, P = 0.012) and structural connectivity (T = 2.708, P = 0.006) outperforming functional synchronization (Fig. 3).

Fig. 3.

Effects of imaging modalities (left) and subnetworks (right) on the prediction performance for the changes of total Beck Scale of Suicidal Ideation scores. DAN, dorsal attention network; DMN, default mode network; FBN, functional brain network; FPN, frontoparietal control network; LN, limbic network; MBN, morphological brain network; SBN, structural brain network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network; *, P < 0.05; **, P < 0.01.

Effects of subnetwork. The SMN was found to outperform the other subnetworks for predicting the percentage changes of total BSSI scores (T = 2.309, P = 0.010) (Fig. 3). However, this difference did not survive correction for multiple comparisons.

Prediction for changes of total HAMD-24 scores

Morphological similarity. The GMV-based morphological similarities between the SMN and DAN (r = 0.491, P = 0.002) and between the SMN and LN (r = 0.597, P < 0.001), the CT-based morphological similarities between the SMN and DAN (r = 0.736, P = 0.001) and between the DAN and LN (r = 0.604, P = 0.001), and the FD-based morphological similarities between the DAN and VAN (r = 0.580, P < 0.001) significantly predicted the percentage changes of total HAMD-24 scores (Fig. 4). Among the GMV-, CT-, and FD-based predictions, the mostly contributed regions were the right median cingulate and paracingulate gyri (Z = 4.063), the left orbital frontal cortex in the LN (Z = 5.119), and the left frontal operculum and insula in the VAN (Z = 2.837), respectively (Fig. 5). No significant predictions were observed for the GI-, or SD-based morphological similarities (P > 0.05, FDR corrected).

Fig. 4.

Prediction performance for the changes of total 24-item Hamilton Depression Rating Scale scores. *, significant linear correlation between real scores and predicted values (P < 0.05, FDR corrected). CT, cortical thickness; DAN, dorsal attention network; DMN, default mode network; FA, fractional anisotropy; FD, fractal dimension; FPN, frontoparietal control network; GI, gyrification index; GMV, gray matter volume; LN, limbic network; rs-fMRI, resting state-functional magnetic resonance imaging; SD, sulcal depth; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.

Fig. 5.

Connectional contributions to the predictions for the changes of total 24-item Hamilton Depression Rating Scale scores. Regional sizes are proportional to total absolute contributions across the linked connectivity. CT, cortical thickness; DAN, dorsal attention network; DMN, default mode network; FA, fractional anisotropy; FD, fractal dimension; FPN, frontoparietal control network; GMV, gray matter volume; LN, limbic network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.

Structural connectivity. The fiber number-based structural connectivity between the FPN and DMN (r = 0.649, P < 0.001), the FA-based structural connectivity between the VN and DMN (r = 0.511, P = 0.002), between the LN and DMN (r = 0.443, P = 0.006), within the VN (r = 0.436, P = 0.007), and within the FPN (r = 0.571, P < 0.001), and the fiber length-based structural connectivity between the VN and DMN (r = 0.703, P < 0.001), between the LN and DMN (r = 0.467, P = 0.004), within the VN (r = 0.494, P = 0.002), and within the FPN (r = 0.449, P = 0.006) significantly predicted the percentage changes of total HAMD-24 scores (Fig. 4). The right precuneus in the FPN (Z = 1.710) contributed the most to the fiber number-based prediction and the left middle and inferior middle temporal gyrus in the DMN contributed the most to both the FA-based (Z = 2.707) and fiber length-based (Z = 2.913) predictions (Fig. 5).

Functional synchronization. The functional synchronization between the VN and FPN significantly predicted the percentage changes of total HAMD-24 scores (r = 0.632, P < 0.001) (Fig. 4). The left calcarine fissure and surrounding cortex in the VN (Z = 4.554) contributed the most to the prediction (Fig. 5).

Effects of imaging modality. Significant effects of imaging modality were observed on the model performance of predicting the percentage changes of total HAMD-24 scores (F = 6.551, P = 0.001), with structural connectivity outperforming morphological similarity (T = 3.615, P < 0.001) (Fig. 6).

Fig. 6.

Effects of imaging modalities (left) and subnetworks (right) on the prediction performance for the changes of total 24-item Hamilton Depression Rating Scale scores. DAN, dorsal attention network; DMN, default mode network; FBN, functional brain network; FPN, frontoparietal control network; LN, limbic network; MBN, morphological brain network; SBN, structural brain network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network; *, P < 0.05; **, P < 0.01.

Effects of subnetwork. The SMN was found to outperform the other subnetworks for predicting the percentage changes of total HAMD-24 scores (T = 2.181, P = 0.014) (Fig. 6). However, this difference did not survive correction for multiple comparisons.

Discussion

This study constructed the predictive models of the treatment response to accelerated NH-rTMS in MDD patients with suicidal ideation based on multimodal large-scale brain networks. We found that the morphological and structural brain networks predicted the percentage changes of both total BSSI and HAMD-24 scores. Additionally, the functional brain networks predicted the percentage changes of total HAMD-24 scores. Further analyses revealed that the structural brain networks showed better prediction performance than the morphological and functional brain networks and the SMN exhibited superior predictive capabilities than the other subnetworks. These findings may help determine personalized treatment plans for MDD patients with suicidal ideation.

We found that the GMV-, CT-, and FD-based morphological brain networks and fiber length-based structural brain networks were able to predict the percentage changes of both total BSSI and HAMD-24 scores and the fiber number- and FA-based structural brain networks and functional brain networks were able to predict the percentage changes of total HAMD-24 scores. These findings suggest that large-scale brain networks may serve as predictive markers for the efficacy of accelerated NH-rTMS treatment in MDD patients with suicidal ideation. Thus, large-scale brain networks may help provide potential clinical guidance for identifying MDD patients who are most likely to benefit from NH-rTMS treatment and thus enhance clinical efficiency by facilitating personalized treatment strategies. Specifically, for a MDD patient with suicidal ideation, we can first collect his or her multimodal MRI data and construct corresponding brain networks. The brain networks can then be fed to the predictive models constructed in this study to predict his or her response to NH-rTMS treatment. Based on the predicted response, we can determine whether the patient is suitable for NH-rTMS treatment. It should be noted that different sets of subnetworks or features contributed to the predictions between the BBSI and HAMD-24 scores. The BSSI measures the severity of suicidal ideation (Beck et al., 1979) and the HAMD-24 quantifies the severity of depressive symptoms (Hamilton, 1967). Therefore, the observed discrepancy in the predictive features between these two measures suggests the distinct neural mechanisms underlying the antidepressant and anti-suicidal effects of accelerated NH-rTMS. In addition, we noted that the GI- and SD-based morphological brain networks could not predict the percentage changes of neither total BBSI nor HAMD-24 scores. Previous studies have demonstrated feature-dependent connectivity patterns and neurobiological correlates for cortical morphological brain networks (Li et al., 2021, 2023). Hence, the failure of the predictions based on the GI- and SD-based morphological brain networks may be due to the insensitivity of brain connectivity information captured by GI and SD to the anti-suicidal and antidepressant mechanisms of accelerated NH-rTMS.

We found that structural networks exhibited better capabilities than morphological and functional networks in predicting both the percentage changes of total BSSI and HAMD-24 scores. Previous studies have shown that DTI-based measures are able to predict rTMS efficacy of antidepressants and potentially identify people who do not respond to rTMS treatment (Chen et al., 2022; Klooster et al., 2020; Ning et al., 2022). In particular, accumulating evidence suggests that TMS pulses travel along white-matter tracts, which influence the propagation of action potentials (Kearney-Ramos et al., 2018; Momi et al., 2021). This may explain why structural connectivity performed well in predicting accelerated NH-rTMS treatment efficacy in MDD patients with suicidal ideation. Based on our findings, we argue for adopting the DTI-guided, network-informed accelerated NH-rTMS to enhance target engagement in future treatment of MDD patients with suicidal ideation. This sounds quite plausible given that structural connectivity is relatively stable on short time scales (e.g., minutes to hours). Compared with structural connectivity, functional synchronization is highly dynamic and state-dependent (Stoffers et al., 2015). This nature may increase interindividual variability and reduce test-retest reliability of functional synchronization. In addition, although functional synchronization is widely accepted to be shaped by structural connectivity (Fotiadis et al., 2024; Z. Wang et al., 2015), there is no one-to-one correspondence between them (Honey et al., 2009). All these factors may explain the weaker predictive power of functional synchronization. To test this speculation, future studies could investigate whether functional synchronization measured during naturalistic fMRI paradigms, which has been previously demonstrated higher test-retest reliability than that estimated during resting state (Wang et al., 2017; Zhang et al., 2022), is associated with improved predictive power for the treatment response to NH-rTMS in MDD patients with suicidal ideation. Finally, as for morphological similarity, both human and animal studies have shown that it is correlated with structural connectivity (Li et al., 2024; Wang et al., 2024). Nevertheless, this correlation is relatively weak, suggesting that morphological similarity actually offers a unique framework for a broader understanding of the brain network organization. Further studies on the similarities and differences between morphological similarity and structural connectivity are needed to help elucidate the weaker predictive power of the former than the latter observed in this study.

In addition to the imaging modality-dependent prediction performance, we further found that the SMN outperformed the other subnetworks in predicting both the percentage changes of total BSSI and HAMD-24 scores. The SMN is involved in basic sensory/motor functions, which are reported to show deficits in MDD (Ahern & Semkovska, 2017). Supporting the findings of functional deficits, various neuroimaging studies have found structural and functional alterations in SMN regions in MDD. For example, an ENIGMA-based mega-analysis reported reduced CT in somatomotor regions in adolescent MDD patients (Schmaal et al., 2017). In line with previous literature, our results provide further evidence for a close relationship between the SMN and pathological mechanisms of MDD. Additionally, the SMN is related to emotional reactivity to stress regulation (Goldin & Gross, 2010; Goldin et al., 2009). Considering that the MDD patients included in this study were accompanied with suicidal ideation, stress-related neural mechanisms might presumably contribute to the SMN-dominated prediction. Finally, it should be noted that MDD is characterized by deficits in reorienting attention to salient events. Therefore, the SMN-dominated prediction may suggest the role of accelerated NH-rTMS treatment in repairing brain connections transmitting information about salience and somatosensation (Kaiser et al., 2015; Korgaonkar et al., 2020). Future longitudinal neuroimaging studies may provide more insight into this speculation.

This study has several limitations. First, the sample size was relatively small, which challenged the generalizability of the results observed in this study. This issue may be exacerbated by the sample heterogeneity given that both adolescent and adult patients were included in this study. Thus, further larger sample studies with more homogeneous patients are needed to validate the robustness and generalizability of our findings. Second, the predictive models constructed in this study were based on single imaging modality. Given the complementary information provided by brain networks derived from different imaging modalities, it is interesting to examine whether integrating the brain connectivity information from different imaging modalities could may improve the predictive accuracy of accelerated NH-rTMS treatment efficacy. Finally, this study predicted the effects of 7 days’ accelerated NH-rTMS treatment in depressed patients with suicidal ideation. Further studies are needed to examine whether brain networks are predictive of longer-term outcomes after accelerated NH-rTMS treatment.

Conclusion

The antidepressant and anti-suicidal effects of accelerated NH-rTMS can be predicted by large-scale brain networks but in an imaging modality- and subnetwork-dependent manner.