Critical role of rhythms in prefrontal transcranial magnetic stimulation for depression: A randomized sham-controlled study.

Repetitive transcranial magnetic stimulation (rTMS) is an alternative treatment for depression, but the neural correlates of the treatment are currently inconclusive, which might be a limit of conventional analytical methods. The present study aimed to investigate the neurophysiological evidence and potential biomarkers for rTMS and intermittent theta burst stimulation (iTBS) treatment. A total of 61 treatment-resistant depression patients were randomly assigned to receive prolonged iTBS (piTBS; N = 19), 10 Hz rTMS (N = 20), or sham stimulation (N = 22). Each participant went through a treatment phase with resting state electroencephalography (EEG) recordings before and after the treatment phase. The aftereffects of stimulation showed that theta-alpha amplitude modulation frequency (fam ) was associated with piTBS_Responder, which involves repetitive bursts delivered in the theta frequency range, whereas alpha carrier frequency (fc ) was related to 10 Hz rTMS, which uses alpha rhythmic stimulation. In addition, theta-alpha amplitude modulation frequency was positively correlated with piTBS antidepressant efficacy, whereas the alpha frequency was not associated with the 10 Hz rTMS clinical outcome. The present study showed that TMS stimulation effects might be lasting, with changes of brain oscillations associated with the delivered frequency. Additionally, theta-alpha amplitude modulation frequency may be as a function of the degree of recovery in TRD with piTBS treatment and also a potential EEG-based predictor of antidepressant efficacy of piTBS in the early treatment stage, that is, first 2 weeks.

INTRODUCTION

Treatment‐resistant depressive (TRD) patients, who have poor medical outcomes, account for approximately one‐third of major depressive disorder (MDD) patients (Rush, 2007; Voineskos, Daskalakis, & Blumberger, 2020). A number of alternatives to pharmacological treatments, including repetitive transcranial magnetic stimulation (rTMS) and an alternative rTMS protocol of intermittent theta burst stimulation (iTBS), have been approved as effective treatment for TRD by the FDA in the United States and investigated in numerous studies (Bakker et al., 2015; Loo, 2008; Loo, McFarquhar, & Mitchell, 2008). Recently, one study with a large number of participants demonstrated that with an FDA‐approved standard protocol, 3,000 pulses of 10 Hz rTMS and 600 pulses of iTBS had equivalent antidepressant efficacy. For these protocols, iTBS worked with just 3 min of stimulation and was more efficient than 37.5 min of 10 Hz rTMS stimulation (Blumberger et al., 2018). In addition, our previous studies have shown that prolonged iTBS (piTBS; 1,800 pulses), which is three times longer than the standard iTBS protocol, as well as 10 Hz rTMS had more effective antidepressant efficacy than a sham condition following a 2‐week treatment session (Li et al., 2014, 2020). In addition, a large double‐blind, sham‐controlled study has demonstrated that the responders who received rTMS showed primary enhancement during the first 2 weeks (O'Reardon et al., 2007), although previous studies used 4–6 weeks rTMS treatment as a standard protocol (Blumberger et al., 2018; van Eijndhoven et al., 2020).

Based on our previously established parameters, the underlying mechanisms of piTBS and rTMS were investigated in several neural imaging studies (Li et al., 2013, 2018). Nevertheless, the electrophysiological neural correlates of treatment effectiveness under the established parameters remain unclear. Furthermore, despite the efficacy of this alternative treatment in TRD, more than half of TRD patients did not respond to the treatment (Blumberger et al., 2018; Li et al., 2014, 2020). As a result, the investigation of the neural correlates and biomarkers of effect for rTMS and iTBS treatment on MDD are a crucial issue in order to enhance the response rate by allowing more efficient selection of treatment protocol (Daskalakis, 2014; Dhami et al., 2021; Kim, Blumberger, Downar, & Daskalakis, 2021).

A growing body of evidence has demonstrated that oscillatory brain activity could interact and synchronize with the externally imposed frequency during (i.e., online) noninvasive brain stimulation, that is, rTMS (for review see Thut et al., 2017). Thut and colleagues showed augmentation of alpha oscillations over the attention‐associated alpha generated in brain regions after entraining alpha activity with rTMS compared to sham conditions in healthy people (Thut et al., 2011). The evidence suggested that endogenous oscillatory neural activities can align with the driving force's rhythm. In addition to the online effects, numerous studies have investigated the longer lasting aftereffects on brain oscillations after the end of the stimulation in both healthy people and those with depression with rTMS treatment (for reviews, see healthy population: Veniero, Vossen, Gross, & Thut, 2015; depression: Noda et al., 2015). Results in depression studies showed that oscillatory activity changed at the stimulation frequency. However, some studies demonstrated broad frequency band modulations. The heterogeneity of depression is likely to be a critical reason for the inconsistent results.

In addition to the uncontrollable innate characteristics of depression, another possible explanation for inconsistent findings is limitations of the analytical methods. These are often based on assumptions of linearity in the data, whereas brain oscillations are both complex and composed of both linear sinusoidal oscillations and nonlinear amplitude‐modulated signals. Although wavelet and Hilbert spectral analysis can resolve the nonstationary signals by adding temporal variation to adaptively determine phase function, and so obtain the instantaneous frequency (Flandrin, 1998; Huang et al., 2009), intermode amplitude modulations are generally overlooked (Huang et al., 2016). For example, Fourier spectral analysis assumes the amplitude is constant and therefore the information of the variation in amplitude modulations is unexploited, meaning conventional data decomposition methods are insufficient to parse the nonlinear information in neural activity (Breakspear, 2017; Buzsáki & Mizuseki, 2014; Nguyen et al., 2019).

In order to investigate amplitude modulation in neural oscillations, Huang et al. (2016) use the Holo‐Hilbert spectral analysis (HHSA) method to provide full informational representations which include carrier frequency (f c), which is the frequency in the conventional frequency analysis, and amplitude modulation frequency (f am) measures. HHSA has been applied in analyzing brain electroencephalography (EEG) signals in visual perception (Juan et al., 2021; Nguyen et al., 2019) and visual working memory (Liang, Tseng, Yeh, Huang, & Juan, 2021) in our previous studies (see Data S1, Supporting Information). Based on this, the analytical method applied in the present study was primarily HHSA, which has not previously been employed in the study of electrophysiological biomarkers in depression.

In sum, in addition to the varied results of brain oscillations induced by rTMS in depression and the limitation of previous analytical methods, less clear evidence for the electrophysiological neural correlates of iTBS effects on brain oscillations in TRD has been provided. Therefore, the present study investigated the neural correlates of the effects of rTMS and iTBS based on the previously established parameters as revealed by resting state brain oscillations of individuals with TRD through applying a novel analytical method, HHSA. This aimed to further investigate potential optimal biomarkers for facilitating therapeutic efficacy in TRD treatment. Through repetitive stimulation, the effects on plasticity of the brain may be lasting after the stimulation. Therefore, the changes in brain oscillations after repetitive modulation using TMS was expected to be associated with the stimulating frequency. Moreover, the amplitude modulation frequency of TBS effect was expected to be possible to investigate by means of HHSA. Specifically, alpha frequency activity was expected to be elicited after giving stimulation at alpha frequency, that is, using 10 Hz rTMS. Theta amplitude modulation frequency was anticipated to be observed after stimulating within the theta range using iTBS.

MATERIALS AND METHODS

Participants

The present study contained part of the data from Li et al. (2020), which included pre‐ and post‐treatment EEG signals. Adult patients, aged from 21 to 70, who were diagnosed with a recurrent major depressive disorder, based on DSM‐IV criteria (American Psychiatric Association, 1994), were recruited. The eligible participants had also failed to respond to at least one antidepressant treatment in the current depressive episode. In addition, participants were required to score at least four on a Clinical Global Impression‐Severity scale (CGI‐S) as well as at least 18 on the 17‐item Hamilton Depression Rating Scale (HDRS‐17) (Hamilton, 1960). Participants with bipolar I or II disorders, who had a history of psychotic disorders, personality disorders, or neurological disorders (e.g., stroke or seizure), with brain implants or cardiac pacemakers, or who were pregnant were excluded. (Regarding the sample size determination, please refer to the Data S1).

The present study was executed according to the Declaration of Helsinki and was approved by the local ethics review committee. Informed consent was provided by each participant. This research was preregistered in the University Hospital Medical Information Network Clinical Trials Registry (Registration number: UMIN000020892), which is a member of the World Health Organization registry network as well as one of the largest clinical trial registries in Asia.

Procedure

The procedure of the experiment consisted of three segments (see Figure 1). First of all, participants went through a diagnostic interview by a medical psychiatrist in order to confirm that participants met the recruitment criteria (the scores of the assessment at Week 0). Participants also underwent brain structural imaging scanning (a T1 structural scan) by magnetic resonance imaging (MRI), which was used for locating the site for brain stimulation (left dorsolateral prefrontal cortex, lDLPFC). Second, a 5‐min resting EEG with eye‐closed was recorded. After that, participants were randomized into three rTMS treatment groups (Group A: prolonged iTBS, piTBS; Group B: 10‐Hz rTMS; Group C: sham) with an equal ratio of participants for each group. For randomization of the treatment assigned to the patients, a computer‐generated random number list was conducted with a block size of 12 by an independent research assistant. Specifically, the random sequence of each block was composed of four sites each in Group A (piTBS), Group B (rTMS) and Group C (sham). The sequence list was concealed on the research assistant's personal computer until the treatment was assigned. All participants then took part in an acute treatment phase with five consecutive days, with one session per day, over 2 weeks for total of 10 sessions. Finally, 5‐min resting EEG with eyes‐closed was recorded, and the assessments, including HDRS‐17 and CGI‐S, were again evaluated (Week 2). The interval between the last treatment session and resting EEG recording was at least 1 day. The clinical trial was a double‐blind study in which neither the patients nor the doctors knew which treatment the patients were receiving and those who executed the experiment were the only individuals who knew the treatment that each participant would receive. Thus, the prospective effect for participants as well as any bias when doctors evaluated the clinical outcomes of patients were avoided, enhancing the reliability of the results or any effects seen. The study took place at the TMS laboratory of Taipei Veterans General Hospital.

FIGURE 1

The procedure of the experiment. The first segment at Week 0, the participants were assessed by CGI‐S and HDRS_17 and underwent brain structural imaging scanning by MRI. Second, at the pre‐treatment stage, a 5‐min resting electroencephalography (EEG) with eye‐closed was recorded. After that, participants were randomized into three rTMS treatment groups (either piTBS, 10 Hz rTMS or sham) and received the treatment with five consecutive days, with one session per day, over 2 weeks for total of 10 sessions. At the post‐treatment stage, a 5‐min resting EEG with eye‐closed was recorded again. Finally, CGI‐S and HDRS_17 were applied to evaluate the clinical symptoms of participants. CGI‐S, Clinical Global Impression‐Severity scale; HDRS_17, 17‐item Hamilton Depression Rating Scale; MRI, magnetic resonance imaging; piTBS, prolonged iTBS

Brain stimulation (TBS/TMS) parameters and site localization

The adopted parameters were based on our previous studies, which showed the antidepressant efficacy in 2‐week treatment sessions (Li et al., 2014, 2016, 2020). TBS and rTMS were delivered by Magstim Rapid2 stimulator (Magstim Co., Ltd., Whitland, UK) connected to a double 70 mm figure‐of‐eight air‐cooled coil. For the TBS, the basic stimulation parameters were three pulses of TMS delivered at 50 Hz and repeated every 200 ms (5 Hz). The prolonged iTBS was adopted based on the standard iTBS protocol (a 2‐s train of TBS repeated every 10 s at an intensity of 80% active motor threshold and a total of 600 pulses) by increasing the number of pulses to 1,800. The rTMS stimulation was a 4‐s train of 10 Hz rTMS repeated every 30s and with a stimulation intensity of 100% resting motor threshold, for a total of 1,600 pulses for one session. For the sham group, a sham coil was used (Magstim Placebo Coil; Magstim Co., Ltd.) to produce similar pulse sounds and mimic the sensation of the pulses while without stimulating the cortex. Participants in Group A received piTBS (1,800 pulses/session × 10 sessions); Group B received 10‐Hz rTMS (1,600 pulses/session × 10 sessions) and Group C received sham treatment (randomly assigned the parameters of piTBS or 10 Hz rTMS, for a total of 10 sessions).

The localization of the left DLPFC was carried out following the methods used in our previous studies (Li et al., 2014, 2020). The left DLPFC was determined as a landmark between the junction of Brodmann area (BA) 9 and 46 according to each participant's brain structural image. For MRI navigation, the structural T1‐weighted images were obtained from a 3.0 GE Discovery 750 whole‐body high‐speed imaging device (GE Healthcare, Milwaukee, WI). Based on participants' T1 image, a brain‐navigation computer and a Polaris infrared tracking system (Brainsight, Rogue Research, Inc., Montreal, QC) were employed to localize the stimulating site at the left DLPFC.

Electrophysiological recording parameters

A Neuroscan amplifier (NuAmps) was used, and a 32‐channel EEG cap (Quik‐Cap) with Ag/AgCl gel was placed on participants' head, positioned according to the international 10‐20 system for recording EEG signals (Channels: FP1, FP2, F7, F3, Fz, F4, F8, FT7, FC3, FCz, FC4, FT8, T3, C3, Cz, C4, T4, TP7, CP3, CPz, CP4, TP8, T5, P3, Pz, P4, T6, O1, Oz, O2, A1, A2). The entire 5‐min eye‐closed resting EEG was digitized at a 1,000 Hz sampling rate without applying online filters. The reference was the average of electrodes at the two sides of the mastoid (A1 and A2). In addition to the electrodes on the scalp, two pairs of bipolar electrodes were mounted in order to detect eye movements with the VEOU and VEOL electrodes placed above and below the left eye, respectively, while the HEOR and HEOL electrodes were positioned adjacent to the canthus of each eye. The impedances of all channels were kept below 5 kΩ.

Efficacy assessments

The symptomatic ratings were taken at baseline (Week 0, W0, before the acute TMS treatment), at the end of first week (W1) of stimulation, after finishing the treatment (W2) as well as at a 12‐week follow‐up (W14) (for the changes of intervention at W14 refer to Data S1). Due to the scores of each participant varying at baseline the efficacy was measured as the percentage change of HDRS‐17 score (% HDRS‐17) from before (W0) and after 2 weeks of treatment (W2) across the three groups. The response rate was when the percentage change of HDRS‐17 score was at least a 50% reduction at W2 compared to the score at baseline, W0, who are called responder.

Data analysis

Behavioral data

The statistical analysis of demographic and clinical data was conducted using PASW Statistics, Version 18.0 (SPSS, Inc., Chicago, IL). Continuous variables such as the scores of assessments across groups were compared using one‐way ANOVA and categorical variables such as gender and response rate among groups were analyzed using Pearson's chi‐square test. Percentage changes of HDRS‐17 at Week 2 for each group and for differences between piTBS, rTMS and sham group were calculated. Bonferroni correction was applied to post hoc analyses. The level of statistical significance was set at p < .05 (two‐sided). The mean differences expressed between groups, denoted as the effect sizes (Cohen's d), were calculated.

Resting EEG‐Holo‐Hilbert spectral analysis

The continuous EEG data were epoched to 10 s for each segment with the first and last 10 s of data excluded to prevent inclusion of unstable data. Afterward, a low pass filter, set at 100 Hz, was applied to eliminate any potential higher frequency contamination from the environment. The baseline was executed to correct the shifting of the signals, which was set at minus infinite to infinite intervals. Afterward, independent component analysis (ICA) was conducted for each participant to eliminate the associated ocular and muscle artifacts. Trials were removed with an artifact rejection criterion of ±200 μV to remove severely drifting and fluctuating signals. The data then was de‐trended and normalized by dividing by the standard deviation of the time serious data. The steps above and statistical analysis were performed using SPM12 for MEG/EEG (Wellcome Department of Cognitive Neurology, London, UK).

After the preprocessing of EEG signals, Holo‐Hilbert spectral analysis was conducted with customized MATLAB (MathWorks) scripts. Two‐layer empirical mode decomposition (EMD) was used for the HHSA. EMD was first introduced by Huang et al. (1998) and is the process of decomposing the broad band signals into several narrow band signals, named intrinsic mode functions (IMFs). The process of EMD involves the repeating sifting process. First of all, the upper maxima and lower minima were determined from the signal. Second, connecting each determined data point by cubic spline, which was used to produce the upper and lower envelope. Third, the mean of the upper and lower envelope was calculated and subtracted from the original signal. The residual signal is the first IMF. The first IMF then becomes an original signal to repeat the aforementioned three steps and afterward obtain the second IMF. This process was repeated until the trend was shown. By this sifting process, the signals would be decomposed to several IMFs from high frequency to low frequency.

Regarding two‐layer processes of HHSA, in the first layer EMD, the EEG signals were decomposed into several intrinsic mode functions (IMFs) with the characteristics of decreasing frequencies from high frequency to low frequency. The second layer IMFs were then obtained by applying EMD to the envelope with absolute values of each first layer IMF. The involved EMD method is masking EMD, which is an improved algorithm proposed by Tsai and colleagues to prevent the mode‐mixing problem (Tsai, Fan, Lin, Huang, & Yeh, 2016). The decomposed IMFs and the distribution of IMFs are shown in Figure 2. Furthermore, the instantaneous carrier frequencies as well as the amplitude of first layer IMFs were extracted by Quadrature and, in addition, the instantaneous amplitude modulation frequencies were obtained from second layer IMFs (Huang et al., 2016). This extracted information was then projected to a space to form the three‐dimensional (amplitude modulation frequency × carrier frequency × time) Holo‐Hilbert spectrum. In the present study, time information was not of concern in the resting EEG data and so, in order to visualize the results, the time dimension of the spectral power was summed to produce a two‐dimensional Holo‐Hilbert spectrum (amplitude modulation frequency × carrier frequency). In this spectrum, the y‐axis represents amplitude modulation frequency and the x‐axis refers to carrier frequency. All the trials were then averaged and the data from the same stimulation group were merged as one dataset. Afterward, the averaged and merged data of each group was logarithmically rescaled to elevate the homogeneity and which could more readily fit to a normal distribution for further statistical analysis.

FIGURE 2

A demonstration of intrinsic mode functions (IMFs). (a) demonstrates 1 s of raw resting EEG data from the Cz electrode within one trial. (b) Shows the mask empirical mode decomposition (EMD) results of one trial from IMF one to IMF ten. Distinguished IMFs refers to the electrophysiological activities in corresponding frequencies. (c) Illustrates the probability of IMF frequency distributions across all participants. The dashed line denotes the partition between two frequency distributions. The cross symbol indicates the peak frequency of each IMF. These were based on the EMD results

For the statistical analyses, the log power changes of each subgroup of pre‐ and post‐treatment, that is, piTBS_Responder, piTBS_NonResponder, rTMS_Responder, rTMS_NonResponder in the Holo‐Hilbert topographic maps were contrasted with the sham group by using independent t tests. A cluster‐based nonparametric permutation (CBnPP) test under p < .05 with 5,000 permutations was used (Groppe, Urbach, & Kutas, 2011; Maris & Oostenveld, 2007). Due to the limitation of the discrepancy of the number of subjects between subgroup of the responder and sham, the level of statistical significance was set at t value >2.048 for piTBS_Responder comparison (degrees of freedom = 28) and t value >2.056 for rTMS_Responder comparison (degrees of freedom = 26) based on the t distribution table. The neighboring distance between two EEG sensors was set at 60 mm. Afterward, post hoc analysis was done by choosing the frequency ranges and brain regions with noticeable power enhancement observed in the Holo‐Hilbert topographic distribution under each active TMS treatment. The power within the chosen range was averaged as a value with arbitrary units. The independent t test was conducted with p‐value setting at p < .05. Then the Pearson's correlation was performed for power changes with the frequency ranges of averaged log power values and the percentage of the HDRS improvement (HDRS_% improvement) in each active treatment group (piTBS, 10 Hz rTMS) to investigate the frequency bands associated with antidepressant efficacy.

RESULTS

Demographic data and antidepressant efficacy of piTBS/rTMS outcomes

All participants in the piTBS (N = 19), 10 Hz rTMS (N = 20), and sham (N = 22) groups completed the procedure within around two and half years (Figure S1). Demographic variables and the scores of clinical assessments at baseline were not different across the three groups (Table 1).

TABLE 1

Demographics and clinical assessments across treatment groups

	piTBS (N = 19)	rTMS (N = 20)	Sham (N = 22)
Age, years	48.74 ± 14.41 a	49.10 ± 14.78	48.50 ± 12.87
Female	14 a	14	16
CGI‐S (BL)	4.32 ± 0.58 a	4.60 ± 0.82	4.41 ± 0.59
HDRS‐17 (BL)	22.53 ± 3.17 a	22.60 ± 3.33	22.59 ± 2.61
% change at W2	−40.85 ± 6.70 b	−30.19 ± 6.26	−14.75 ± 2.94
Comparing with sham, p value	<0.01**	0.53
Responders at W2	8 (42%) a	6 (30%)	0 (0%)**

Abbreviations: BL, baseline; CGI‐S, the Clinical Global Impressions Scale‐Severity of Illness; HDRS‐17, 17‐item Hamilton Depression Rating Scale; piTBS, prolonged intermittent theta burst stimulation; rTMS, repetitive transcranial magnetic stimulation; W, week.

Values are mean ± SD, N or N%.

Values are mean ± SE.

**p < .01.

The antidepressant efficacy results demonstrated that the percentage changes of HDRS‐17 at W2 were significantly different between the three cohorts of treatment, F (2,58) = 6.00, p < .01. The post hoc test showed that the antidepressant efficacy of piTBS was significantly larger than the sham group effect (M = −40.85%, SE = 6.70 vs. M = −14.75%, SE = 2.94; p < .01; Cohen's d = 1.20). Although rTMS showed larger antidepressant efficacy in comparison with sham, the effect did not reach statistical significance (M = −30.19%, SE = 6.26 vs. M = −14.75%, SE = 2.94; p = .13; Cohen's d = 0.73). The percentage changes of the depression score at W2 did not show a difference between piTBS and rTMS (M = −40.85%, SE = 6.70 vs. M = −30.19%, SE = 6.26; p = .53; Cohen's d = 0.38) (Figure 3). Additionally, the response rate was significantly different between the three groups, χ 2 (2, N = 61) = 11.06, p < .01. piTBS clinical outcomes showed a 42% (8 out of 19 patients) response rate and rTMS showed a 30% (6 out of 20 patients) response rate. No responders at Week 2 were observed in the sham group.

FIGURE 3

The antidepressant efficacy of TMS treatment outcome in each group. The percentage changes of HDRS‐17 at week two were significantly larger in piTBS compared to the sham group (p < .01). However, the antidepressant efficacy in rTMS did not significantly differ from piTBS or the sham group. The results indicate piTBS had antidepressant efficacy on TRD, while rTMS did not have such a notable treatment effect in the sub‐dataset applied in the present study, which was from Li et al. (2020). **p < .01. W = week; Error bars are standard errors (SE)

Generally, all participants tolerated the treatment well. Some reported temporary headaches (piTBS, N = 2; 10 Hz rTMS, N = 3; sham, N = 1, χ 2 (2, N = 61) = 1.31, p = .52), temporary dizziness (piTBS, N = 3; 10 Hz rTMS, N = 4; sham, N = 3, χ 2 (2, N = 61) = 0.32, p = .85), and one patient who received 10 Hz rTMS treatment reported tinnitus. No patients experienced seizure or mania.

EEG oscillation results

EEG data from all participants were subject to analysis. To investigate the neural oscillation changes related to treatment efficacy, the changes of power, which is the contrast of pre‐ and post‐treatment resting EEG, in active treatment groups involving piTBS_Responder, piTBS_NonResponder, rTMS_Responder, and rTMS_NonResponder were compared to the sham group, respectively, using independent t tests. From the Holo‐Hilbert topographic distribution of piTBS_Responder in comparison to sham, increased power was observed within the theta‐alpha amplitude modulation range (3–11.8 f am), which modulated the theta‐alpha carrier frequency (3–11.8 f c) over frontal regions and occipital regions (Figure 4a). However, such power changes were not seen for piTBS_NonResponder compared to sham (Figure 4b). From the post hoc analysis, based on the topographic maps, the averaged log power within theta‐alpha amplitude modulation range (3.67–11.31 f am|8.72–11.31 f c) over frontal regions (F3, Fz, F4, FC3, FCz, FC4) was significantly increased in piTBS_Responder relative to sham, t(28) = 2.59, p < .05 (Figure 4c). Furthermore, the changing power of theta‐alpha amplitude modulation frequency in piTBS was found to be positively correlated with the HDRS% improvement, r (19) = .51, p < .05 (Figure 4d). Although the power difference between piTBS_Responder and sham was also seen in occipital regions, the average log power within the same frequency range (3.67–11.31 f am|8.72–11.31 f c) over occipital regions (O1, Oz, O2) did not show a significant correlation with symptom improvement, r (19) = .44, p = .06. The same frequency ranges and brain regions as piTBS_Responder were selected to compare the power difference between rTMS_Responder and sham. No significant difference was found, t(26) = 1.09, p = .29.

FIGURE 4

The subgroup of the piTBS treatment effect compared to sham group in TRD and the following correlation with clinical outcome. (a) Shows the Holo‐Hilbert topographic distribution of different power in piTBS_Responder in comparison to the sham group by cluster permutation independent t test. The x‐axis indicates the carrier frequency and the y‐axis indicates the amplitude modulation (AM) frequency. The color bar represents the t value. The positive t value, represented by the regions of redder color, denotes an increment in power in piTBS_Responder relative to sham group. The negative t value, represented by areas of blue, denotes a decrement in power in piTBS_Responder. The log power spectrum shows the trend of increased power at the 3–11.8 Hz modulated alpha carrier frequency, involving the range of the theta to alpha rhythm that corresponded to the theta burst frequency of the stimulation. Responder denotes that the percentage changes of HDRS‐17 at W2 were larger than 50%. (b) Shows the Holo‐Hilbert topographic map of power change in piTBS_NonResponder in comparison to the sham group by cluster permutation independent t test. The pattern of the distribution was different from piTBS_Responder. No theta‐alpha amplitude modulation has been found in piTBS_NonResponder compared to the sham group. (c) Shows the post hoc analysis that the changes of log power of the theta‐alpha amplitude modulation frequency within the AM range from 3.67 to 11.31 Hz and carrier frequency range from 8.72 to 11.31 Hz over frontal regions in piTBS_Responder compared to sham. The changes of log power of the theta‐alpha amplitude modulation frequency significantly increased in piTBS_Responder, as shown by the independent t test. *p < .05. The error bars indicate the standard error (SE). (d) The Pearson's correlation was calculated to test the relationship between the improvement of the depressive symptoms and the difference of theta‐alpha amplitude modulation frequency log power at frontal regions between pre‐ and post‐piTBS. The results showed a significant positive correlation (p < .05), indicating that the entrainment of theta‐alpha amplitude modulation frequency could be a function of the degree of recovery in piTBS treatment and the biomarker for effective TRD treatment

With respect to EEG power changes in rTMS_Responder and NonResponder compared to sham, the increment of carrier alpha frequencies within the 6.4–11.8 Hz range were both observed across a range of the amplitude modulation frequency bands (0.5–3 f am) (Figure 5a,b). In view of the topographic distribution, the post hoc analysis showed the averaged log power within alpha carrier frequency range (0.5–3.36 f am|8.72–12.34 f c) over frontal‐central regions (F3, Fz, F4, FC3, FCz, FC4, C3, Cz, C4) was significantly increased in rTMS_Responder and rTMS_NonResponder relative to sham, t(26) = 2.07, p < .05 and t(34) = 3.38, p < .05, respectively (Figure 5c,d). Nevertheless, the change in power in the alpha carrier frequency in rTMS did not correlate with the HDRS% improvement, r (20) = .06, p = .81 (Figure 5e). The same frequency ranges and brain regions as rTMS_Responder were selected to compare the power difference between iTBS_Responder and sham. No significant difference was found, t(28) = 1.40, p = .17.

FIGURE 5

The subgroup of the 10 Hz rTMS treatment effect in comparison to the sham group in TRD and the following correlation with clinical outcome. (a) Shows the Holo‐Hilbert topographic distribution of different power in rTMS_Responder by cluster permutation independent t test. The log power spectrum demonstrated an increment in power in the 6.4–11.8 Hz carrier frequencies, which is within the range of the alpha rhythm and corresponded to the 10 Hz TMS stimulation. (b) Shows the Holo‐Hilbert topographic map of power change in rTMS_NonResponder in comparison to sham group by cluster permutation independent t test. The pattern of the distribution was similar to rTMS_Responder. The increment in power was shown in the 6.4–22.6 Hz carrier frequencies. (c, d) Show the post hoc analysis that the changes of log power of the alpha carrier frequency within the AM range from 0.5 to 3.36 Hz and carrier frequency range from 8.72 to 12.34 Hz over frontal‐central regions in rTMS_Responder and rTMS_NonResponder compared to sham, respectively. The changes of log power of the alpha carrier frequency significantly increased in rTMS_Responder and rTMS_NonResponder, as shown by the independent t test. *p < .05. The error bars indicate the standard error (SE). (e) The Pearson's correlation was conducted for log power of alpha carrier frequency at frontal‐central regions and depressive symptom improvement. No significant correlation between the percentage improvement of HDRS and the difference of log power of alpha carrier frequency pre‐ and post‐rTMS was observed. This suggests that alpha carrier frequency may not be a function of degree of recovery in rTMS treatment

DISCUSSION

The present study evaluated specific neurophysiological biomarkers associated with piTBS treatment for TRD and was achieved by employing a novel analytical method, Holo‐Hilbert spectral analysis. This method particularly allows assessment of amplitude modulation, which cannot be measured by traditional EEG spectrum analysis. The augmentation of power changes at the theta‐alpha amplitude modulation frequency and alpha carrier frequency were investigated after delivering TMS in comparison to sham in the related frequencies, that is, piTBS and rTMS, respectively. Theta‐alpha amplitude modulation frequency was shown to be associated with antidepressant efficacy in patients receiving piTBS treatment, while this was not the case for rTMS and the alpha frequency.

Although the clinical outcomes from the subdataset of Li et al. (2020) did not show a larger antidepressant effect in terms of percentage changes of HDRS for rTMS (30.19%) in comparison to the sham (14.75%) which may be due to the smaller effect size, the full dataset from Li et al. (2020) has shown that piTBS and rTMS have a similar antidepressant effect which was significantly larger than sham. This finding was consistent with the results of our previous studies (Li et al., 2021) and similar to Blumberger et al. (2018). In spite of the comparable treatment effect based on the two stimulating protocols, the underlying antidepressant physiological function of rTMS and piTBS might be different and might be represented by distinct patterns of changes in brain oscillations. In addition, the frequency changes, which was correlated with the symptoms improvement, that is, theta‐alpha amplitude modulation, may be as a function of the degree of recovery. However, this effect might be specific to piTBS, not for rTMS due to the different processes in the brain.

Alpha carrier frequency and theta‐alpha amplitude modulation frequency changes were observed in the Holo‐Hilbert spectrum after the stimulation involving the corresponding frequencies, that is, 10 Hz rTMS and piTBS. The 10 Hz rTMS contains alpha frequency which induced a neural entrainment aftereffect in the brain mainly on the alpha carrier frequency without specific amplitude modulation. The possible reasons for the antidepressant efficacy by 10 Hz rTMS may be as follows. Previous studies have shown that the alpha rhythm in the cortex was strongly affected by a thalamic alpha generator (Buzsaki & Draguhn, 2004; Hughes & Crunelli, 2005). Nevertheless, MDD has been envisaged as thalamocortical dysrhythmia (Schulman et al., 2011) and impaired prefronto‐thalamic functional connectivity has been demonstrated in TRD (Li et al., 2013) with findings showing that the alpha oscillation was less dominated by the thalamus in MDD. After prefrontal rTMS treatment, the imposed alpha activity may exert top‐down control in resetting the cortical, thalamo‐cortical alpha oscillators (Leuchter, Cook, Jin, & Phillips, 2013) as well as recovering the prefronto‐thalamic connectivity to modulate the aberrant activities in the brain areas which are associated with mood regulations, such as anterior cingulate cortex (ACC) and anterior prefrontal cortex (Sadaghiani et al., 2010), and consequently mitigate depressive symptoms. Li et al. (2013) illustrated an increased relationship between EEG alpha activity and thalamic metabolism after successful rTMS treatment. Another perspective is based on the inhibition hypothesis of alpha oscillations (Klimesch, Sauseng, & Hanslmayr, 2007), which views elevating the synchrony of alpha oscillations as inhibiting the overactivity in such critical circuits in MDD involving the prefrontal cortex and the ACC (Murray, Wise, & Drevets, 2011; Sellmeijer et al., 2018; Yamamura et al., 2016), and consequently improving symptoms. However, the augmentation of alpha power was not correlated with clinical outcomes in the present study, which may indicate that alpha oscillation changes might not be a critical function of depressed symptom recovering.

Regarding the piTBS antidepressant effect, we found that the log power of theta‐alpha amplitude modulation frequency increased significantly after piTBS, especially in Responders, and the theta‐alpha amplitude modulation frequency positively correlated with antidepressant effects (Figure 4). HHSA clearly reflected the entrainment aftereffect of theta burst stimulation. The cross‐frequency interaction showed theta‐alpha amplitude modulation modulated alpha carrier frequency. Our previous imaging study showed that theta burst stimulation primarily modulated the fronto‐cingulate circuit, including ACC and medial PFC (i.e., two theta‐prominent brain regions), in patients with MDD (Li et al., 2018). Previous studies have also suggested that resting‐state ACC theta activity before antidepressant treatment could predict treatment responses in patients with MDD (Pizzagalli, 2011; Pizzagalli et al., 2018). In addition, higher ACC theta activities were reported in responders than nonresponders (Korb, Hunter, Cook, & Leuchter, 2009; Pizzagalli et al., 2001). Furthermore, frontal EEG theta from the scalp could be augmented by a cognitive task and was found to be associated with more ACC glucose uptake, which could also predict better antidepressant outcomes (Li et al., 2016). According to Buzsaki (2006), high frequency neuronal activities are prone to decrease rapidly and are more local. However, lower frequency oscillations could influence more neurons, be more long‐lasting and act longer distance. The pattern of the present results showed the nest of fast (alpha) oscillations into slower (theta‐alpha) amplitude modulations may illustrate that the amplitude modulation of theta oscillations might propagate to the main circuits associated with MDD, result in long‐range network synchronization. The patients who had significant enhancement of theta‐alpha amplitude modulation frequency after piTBS showed more marked pathological improvement, which might indicate that the theta‐alpha amplitude modulation frequency changes could serve as a function of the degree of recovery by piTBS. Furthermore, this indicates that the interactions of neuronal oscillations across multiple frequencies between theta‐alpha amplitude modulation frequency and rhythms could be a potentially pivotal EEG predictive biomarker for the antidepressant effect for monitoring the antidepressant effect of piTBS in the early treatment session, that is, first 2 weeks. Previous studies generally applied 4–6 weeks rTMS treatment as a standard protocol (Blumberger et al., 2018; van Eijndhoven et al., 2020). From the EEG patterns of theta‐alpha amplitude modulation frequency changes after 2 week's treatment in the present study, the NonResponder did not show the entrainment aftereffects by piTBS, which is clearly dissociated from Responder. This may provide the insight for doctors for the decision that whether keep using the present protocol of stimulation or change to other treatments, which could save time for the patients and make the treatment more efficient. The cross‐frequency oscillatory interactions between theta‐alpha amplitude modulation frequency and alpha carrier frequency were demonstrated by applying high‐dimensional nonlinear analysis in the piTBS group. In addition to the alpha oscillatory effect in the thalamo‐cortical circuit, theta rhythm may influence the core mechanisms relating to the fronto‐cingulate connectivity in MDD. This might be the reason that piTBS has more efficient antidepressant efficacy than rTMS. Additionally, the modulation by iTBS may be at the neurochemical level, that is, the balance of the GABAergic and glutamatergic systems, which is a critical mechanism in depression (Iwabuchi et al., 2017; Li et al., 2019). Iwabuchi et al. (2017) provided evidence that the ratio of GABA/Glx was reduced and associated with the increased effective connectivity of the DLPFC‐right anterior insula after iTBS compared to sham control in healthy people. This may shed light on the possible mechanism of iTBS in pathophysiology of depression.

The measurements of EEG and TMS were not simultaneous in the present study, which means it may not provide direct evidence of an entrainment effect. The entrainment effect refers to a unidirectional process in which external rhythms may reset the self‐sustained oscillators to oscillations the same as the external frequency and furthermore altering the core brain network activity in associated sensory and cognitive processes (Lakatos, Gross, & Thut, 2019; Thut, Schyns, & Gross, 2011). The concurrent recording of EEG with TMS delivery may be beneficially employed in future studies. However, the worry is that the entraining time may be too short to gain the entrainment effect in the treatment of depression in one session. Additionally, the number of responders in one treatment condition was limited. However, this did not influence our main findings in that the critical biomarkers of brain oscillations at the pretreatment state could differentiate responders and nonresponders. This does, however, require further large‐scale, sham‐controlled studies. In order to increase the response rate, the entrained frequency of theta‐alpha modulation may be tuned to individual rhythms, which may consequently improve the progress of precision medication in the future.

CONCLUSIONS

The present study, utilizing HHSA to study the nonlinear cross‐frequency modulation, provided evidence of the distinct antidepressant neural correlates of 10 Hz rTMS and piTBS for TRD. The findings suggested that the TMS treatment, whatever protocol was applied, may provide lasting benefits, with the possible entrainment aftereffects seemingly revealed in resulting brain oscillation changes. Specifically, while the enhancement of alpha oscillations was associated with 10 Hz rTMS, theta‐alpha amplitude modulation frequency was found in piTBS_Responders. Furthermore, our findings imply that theta‐alpha amplitude modulation frequency changes may indicate a function of response to piTBS treatment, not for rTMS, which may also indicate a possible predictive biomarker for treatment efficacy of piTBS in the early treatment stage, that is, first 2 weeks.

CONFLICT OF INTEREST

The authors declare no potential conflict of interest.

AUTHOR CONTRIBUTIONS

Chi‐Hung Juan and Cheng‐Ta Li designed the experiments. Yi‐Chun Tsai and Cheng‐Ta Li collected the data. Yi‐Chun Tsai, Chong‐Chih Tsai, Norden E. Huang, and Wei‐Kuang Liang analyzed the data. All authors wrote the manuscript.

Appendix S1: Supporting Information (We have provided APA reference style for new version of supporting information as in the attachment)

Click here for additional data file.