A systematic review and NMA were performed according to the guidelines of the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) statement and its extension (18, 19). Various published literature reports that included RCTs were searched using the PubMed, Cochrane Central Register of Controlled Trials, and EMBASE databases from their inception to June 20, 2022. Different frequencies, locations of the coil, intensities (the percentage of the resting motor threshold; RMT), and treatment durations of rTMS to treat adults with TRD were considered while searching the literature. When screening the trials, careful measures were taken to search only selected published literature by applying the key medical subject headings or a combination of their text words like “transcranial magnetic stimulation” and “depressive disorder, treatment-resistant” (Supplementary Table 1). Full-text articles, both peer-reviewed and non-peer-reviewed in English, were retrieved and evaluated by two independent reviewers (JL and HL) who performed the search; any disagreements were resolved through discussion with a third reviewer (LC).

Randomized controlled trials (RCTs) of different parameter selections of rTMS used to treat TRD were included, and the following study characteristics were recorded: (1) Participant: Participants were limited to adults (over 18 years; both females and males) with a primary diagnosis of TRD including unipolar or bipolar depressive episodes and without secondary mood disorders, according to standard diagnostic criteria (DSM-IV, DSM-IV-TR, and ICD-10) (20). TRD was defined for MDD as failing to respond to at least two adequate doses and durations of antidepressant medications (4). (2) Intervention: Trials that have conducted at least two of the following interventions were included: Usually, in clinical practice, two significant modalities of rTMS are applied, i.e., low frequency and high frequency. However, both have specific functional activities; for example, rTMS with a high frequency (HF) and a low frequency (LF) are thought to impart a stimulating and inhibitory effect on the cerebral cortex, respectively. RMT > 100% was defined as sup-threshold rTMS (sup-rTMS), while RMT ≤ 100% defined as sub-threshold rTMS (sub-rTMS). Low frequency is graded as <5 Hz rTMS over the right dorsolateral prefrontal cortex (DLPFC) target at a delivered RMT > 100% (LFR-sup-rTMS). High frequency is graded as ≥5 Hz rTMS over the left DLPFC at a delivered RMT ≤ 100% (HFL-sub-rTMS) (21, 22). HFR-sub-rTMS, HFL-sup-rTMS, LFL-sup-rTMS, LFL-sub-rTMS, and LF stimulation were delivered to the right, followed by HF stimulation over the left DLPFC at the resting motor sup-threshold (LFR-HFL-sup-rTMS), HF over the left followed by LF over the right DLPFC at the resting motor sup-threshold (HFL-LFR-sup-rTMS), LF over the left followed by LF over the right DLPFC at the resting motor sup-threshold (LFL-LFR-sup-rTMS), and the sham control (Figure 1). (3) Comparison: The control group was restricted to another active treatment or the sham control. (4) Outcome: From all included studies, the main focus was to check the response rate of the treatment. Also, evaluation of depressive symptoms in TRD patients was performed with either the Hamilton Depression Rating Scale, 17 items, 24 items, or 28 items (HAMD-17, HAMD-24, and HAMD-28) or the Montgomery Asberg Depression Rating Scale (MARDS). (5) Study design: RCTs and randomized crossover trials were included, but only the first-phase data from crossover trials were obtained (23). Simultaneously, non-randomized controlled trials, quasi-randomized trials, and incomplete trials with only up to 20% of data or enrolled participants with secondary mood disorders like post-stroke depression or psychotic depression with Parkinson’s disease and other concomitant severe medical illnesses were excluded (20). (6) Papers published in a language other than English were excluded (14). The first and third authors (JL and HL) independently selected the studies, and all authors (JL, HL, and LC) discussed and resolved disagreements.

Data were initially gathered by the first and third authors (JL and HL) independently; later, all authors (JL, HL, and LC) scrutinized, discussed, and resolved any disagreements and finalized the studies. The following predefined hierarchical characteristics, which are the essential factors according to the study, are presented in Supplementary Table 2: (1) Sources: Careful screening was done on the gathered data to report the details like the first author, publication year, and so on; (2) Details: Details regarding the ethnicity of the population were studied along with a clear period or duration of the study. Furthermore, demographic and clinical details like the sample size, mean age, and medication details, whether the patients are taking antidepressant drugs or not, intervention measures like rTMS frequency, location of the coil, intension, number of pulses per session, intervention duration, and depression scale were documented carefully; (3) Results: In-depth observation and recording of the minute details regarding the response rate of the treatment, all-cause discontinuation, observing the remission rate by any chance, and the baseline and endpoint scores of the depression scale were graded.

The response and remission rates based on the primary outcome scale were provided. To evaluate the data efficacy and acceptability, the corresponding response rate and discontinuation rate were selected as the primary outcome. Secondary outcomes were the remission rate and the endpoint depression score. The patient’s treatment was considered valid if the response rate improvement was ≥50% from the baseline score of the study depression scale, and remission was considered if the remission rate was found to be ≤7, ≤8, or ≤10 on the HAMD-17, HAMD-21, or MADRS, respectively. The discontinuation rate is the percentage of patients who withdrew from the RCT for any reason.

The risk of bias was assessed according to the Cochrane Handbook for Systematic Reviews of Interventions (version 5.1.0) and the certainty of evidence. The details for judging the risk of bias were as follows: The investigators clearly described random components in the sequence generation process based on the score criteria of low risk, high risk, and unclear risk. Similarly, other criteria for judging the risk of bias included allocation concealment, blinding of participants and personnel, blinding of outcome assessment, incomplete outcome data, selective reporting, and other bias score criteria based on their low, high, and unclear risks. Two authors (JL and HL) independently evaluated the risk of bias and discussed reaching an agreement.

STATA version 15.1 was used to perform the pairwise analysis and NMA, while Revman 5.3.3 was used to perform the risk of bias assessment (24). For all data, the appropriate statistical analysis was applied wherever necessary. The relative risk (RR) for dichotomous outcome measures and standardized mean differences for continuous outcomes, combined with a 95% confidence interval (CI), were estimated. In certain cases, if the number of RCTs included was ≥2, pairwise meta-analyses were conducted through a random-effects model to estimate the direct effects of the interventions. Furthermore, statistical heterogeneity was assessed using I² statistics (25–27). If the estimated I² values were <25%, between 25 and 50%, or >50%, the heterogeneity was considered low, moderate, or high, respectively. If I² was ≥50%, subgroup analysis was performed to explore the source of the heterogeneity (27, 28).

Subsequently, NMA was carried out to analyze the mixed effects of direct and indirect comparison based on the frequentist framework using random-effects models (24–26). The network plots were generated for each outcome (i.e., the response rate, remission rate, discontinuation rate, and endpoint depression score) to clarify the direct or indirect comparison of the included studies. Next, the overall rankings for the efficacy and acceptability of each treatment were estimated using the surface under the cumulative ranking curve (SUCRA), which ranged from 0 to 1. The closer the SUCRA value was to 1, the better the treatment effect (24, 29).

A common heterogeneity parameter (similarity, transitivity, and consistency) was assumed for all treatment comparisons. Similarity and transitivity require combining studies (e.g., A vs. C and B vs. C studies) similar to their effect modifiers (study type and population baselines, such as sample size, mean age, and sex) (30). To illustrate similarity, the RCTs of those groups having no significant differences in the population baseline data were included using eligible trials (16). Consistency was evaluated in the loop-specific approach using the inconsistency factor (IF) to compare the direct and indirect treatment effects (14, 31). In addition, the global Wald and the node splitting approach were used to evaluate the heterogeneity of the network (30–32). A P < 0.05 for the global Wald or the node splitting test indicated the presence of a statistically significant inconsistency in the network model.

Furthermore, subgroup analyses were conducted to explore the difference in the primary outcomes based on different intervention durations between the two groups (divided into the <4 weeks group and ≥4 weeks group). Also, a funnel plot was made to estimate small-study effects, including comparisons of active intervention and a sham and between active intervention and another active intervention (33). Finally, sensitivity analyses, in which studies with too short (one week) of a treatment duration were excluded, were performed to assess the robustness and reliability of the results.

A pairwise meta-analysis was conducted for eight groups of inter-comparisons to compare the response rates. First, the sham control was compared with each of the following groups: sham vs. HFL-sup-rTMS (11 RCTs, 200 participants), sham vs. HFL-sub-rTMS (12 RCTs, 213 participants), and sham vs. LFR-HFL-sup-rTMS (10 RCTs, 118 participants). At a later stage of analysis, LFR-sup-rTMS was compared with HFL-sup-rTMS (2 RCTs, 80 participants) and LFR-HFL-sup-rTMS (2 RCTs, 91 participants). Finally, HFL-sup-rTMS (4 RCTs, 80 participants) was compared with LFR-HFL-sup-rTMS, and HFL-sub-rTMS (2 RCTs, 16 participants) was compared with LFL-sub-rTMS.

More efficacious rTMS interventions than the sham control included HFL-sup-rTMS (RR = 4.35, 95% CI: 2.41–7.87, 12 RCTs, 213 participants), HFL-sub-rTMS (RR = 2.24, 95% CI: 1.03–4.87, 12 RCTs, 222 participants), and LFR-HFL-sup-rTMS (RR = 2.36, 95% CI: 1.26–4.40, 10 RCTs, 118 participants). There was no significant difference between the two rTMS groups.

All comparisons were shown to have a low-to-moderate heterogeneity (I² < 50%), except for HFL sub-rTMS (I² = 68.8%) vs. the sham control. Hence, sensitivity analyses were performed by excluding one of the studies from the HFL-sub-rTMS group vs. the sham control (12 RCTs, 222 participants) at the same time, and the impact of eliminating each of the studies on the overall results between HFL-sub-rTMS and the sham control was estimated.

After the study by Theleritis et al. was excluded, the I² value of HFL-sub-rTMS vs. sham was shown to be 46.9%, with an RR value of 1.72 and a 95% CI of 0.93–3.17 (39). The full text was reviewed, and the risk of bias for the study was reassessed. This study was graded as having an unclear risk of bias. Hence, it was excluded. The detailed information of the sensitivity analyses is reported in Supplementary Table 12, and the direct comparison results, remission rates, endpoint depression scores, and all-cause discontinuation are listed in Supplementary Appendix 4.

The two funnel plots for efficacy and acceptability of treatment are depicted in Supplementary Figure 9. Overall, the funnel plots did not suggest any substantial asymmetry; thus, the presence of a publication bias was not suspected.

Our study has some limitations that must be addressed. First, there was uncertainty in the estimations of the treatment outcomes. A few direct comparisons between two active interventions were included, but 15 of 35 did not report the remission rate. In addition, through 37 RCTs (2 RCTs were not analyzed due to a high risk of bias), 15 of 37 interventions were combined with more than two treatment arms, including 2120 patients.

Second, small-study effects possibly introduced bias as 65% of the trials included fewer than 50 participants. Also, because depressive symptoms in TRD patients were evaluated through either HAMD or MARDS, some methodological heterogeneity might exist. Therefore, the HAMD score followed by the MADRS score was used to calculate the response and remission rates; four studies using only MADRS as the assessment measure were included. Besides, approximately 30% of the studies had an unclear risk of bias, and two studies with a high risk of bias were removed from the analysis. Most of the trials presented an unclear randomization and assignment bias risk, leading to uncertainty for estimates of treatment efficacy.

Additionally, bipolar disorder is a potential cause of TRD (14). Of the included studies, only one trial included patients with bipolar TRD. Thus, subgroup analysis was not conducted to compare the efficacy of rTMS treatment between bipolar and unipolar depression. More trials regarding rTMS treatment of bipolar depression are warranted to gather more evidence. Also, we only used the all-cause discontinuation to assess the acceptability in the present study. However, the higher acceptability may result from efficacy instead of tolerability, leading to an overestimation of tolerability. Therefore, dropout due to the side effects needs to be investigated in future studies.

Another important shortcoming is the failure to assess the degree of cognitive improvement. Cognitive impairment is an important clinical feature of TRD and indicates a poor response to medication. However, only a few of the 37 RCTs conducted pre- and post-treatment cognitive assessments, resulting in insufficient data to allow meta-analysis. Moreover, the treatment duration of two RCTs was only one week. Therefore, sensitivity analysis was performed to assess the robustness and reliability of the results by removing these two trials. The sensitivity analysis results were not affected by these two trials (Supplementary Table 13).

Finally, although we used SUCRA to estimate the ranking order of the comparative effectiveness, the results need to be interpreted with caution because this method may only provide supportive but not conclusive evidence for treatment options (46, 47).