Functional Connectivity Mapping for rTMS Target Selection in Depression

Repetitive transcranial magnetic stimulation (rTMS) to the left dorsolateral prefrontal cortex (DLPFC), a treatment approved by the U.S. Food and Drug Administration for major depressive disorder, induces antidepressant responses in some but not all individuals (1). To optimize treatment, approaches guided by functional MRI (fMRI) have been developed that select target sites based on their functional connectivity (FC) properties (2–10). These methods are increasingly being implemented in experimental rTMS research protocols.

The most common fMRI-guided targeting approach builds on studies associating negative functional connectivity of the stimulation site with the subgenual anterior cingulate cortex (sgACC-StimFC) with better treatment outcomes (11–15). However, strong effects in pioneering studies (3, 5, 6, 8, 16) have been intermixed with negative reports (4, 9, 17, 18). The median sample size for these studies was 25 subjects.

Precise effect size estimates are specifically important for novel biomarkers being used in clinical applications. This report extensively probes possible explanations for mixed associations between sgACC-StimFC and clinical outcomes. First, by examining sgACC connectivity in single subjects with large amounts of fMRI data (“precision functional mapping” data sets), we show that existing methods (2) compensate adequately for noisy sgACC signals. Second, by examining the association between sgACC-StimFC and treatment outcome in the largest data set to date—295 patients with major depressive disorder who received either 10 Hz repetitive transcranial magnetic stimulation (rTMS) or intermittent theta-burst stimulation (iTBS) (19)—we find the model used to estimate the stimulated area in the DLPFC to be a relevant source of variability. Third, by subsampling this large population, which showed a weak overall effect, we demonstrate that sampling variability alone can account for the variable efficacy reported in the literature. Finally, by examining the association between signal quality and treatment prediction, we show, unexpectedly, that a subsample of the population with pronounced global signal variation drives the relationship of sgACC-StimFC to clinical outcome; this subsample was characterized by an especially high variance in global fMRI signals that occur in association with a known pattern of breathing that is prominent in some patients with depression (20).

Collectively, these results provide strong evidence for a true association between FC of the stimulation site with the sgACC and clinical outcome but call into question the utility of current FC-based targeting protocols that rely on single-echo fMRI data and isolated sgACC FC measures, because our results imply that only about 3% of clinical outcome variability is modifiable by this approach.

Methods

Details on all methods and procedures can be found in the online supplement.

Validation of the Weight-Map Approach to Imputing sgACC Signals

Two densely sampled resting-state fMRI data sets were used to validate the weight-map approach (2): a single-echo fMRI data set (named the Midnight Scan Club [MSC] data set in the original publication) of 10 subjects (five female) each with 5 hours of scan time (21) and a multiecho fMRI data set (named the multiecho [ME] data set in the original publication) with five subjects (all male) each with up to 15 hours of scan time (22). See Figure S1 in the online supplement for a graphical explanation of the weight-map method.

Modeling of Stimulation Sites in the DLPFC

Estimating sgACC-StimFC requires a model of the cortical area within the DLPFC that was stimulated in an rTMS session. This is not a straightforward process, as parameters related to stimulation itself and an individual’s anatomy impact the shape of the effective electric field (E-field), which is often multifocal and does not evenly surround the target coordinate. Two competing methods were used to model where rTMS stimulation affected tissue in the DLPFC: E-field modeling performed with the MagVenture B70 coil file in SimNIBS (23) and modeling with the generic “weighted-cone” approximation of the E-field as a 12-mm distance-weighted hemisphere (2–6, 8, 9). For both modeling approaches, sgACC-StimFC was ultimately derived by averaging all of the sgACC FC features within the cortical surface area covered by the relevant E-field or weighted cone after applying thresholds for E-field strength or diameter, respectively.

Modeling Clinical Outcome in THREE-D Patients as a Function of sgACC-StimFC

Sample and procedure.

Details on the THREE-D study (examining the effectiveness of theta burst versus high-frequency rTMS in patients with depression), including the population and TMS treatment protocols, have been published elsewhere (19). Briefly, in the THREE-D study, 414 patients with major depressive disorder were randomly allocated in a 1:1 ratio to receive 20 treatments of either 10 Hz rTMS or iTBS to the left DLPFC. Targeted treatment was based on back-transformation of a predefined stereotactic coordinate (x=−38, y=44, z=26) in Montreal Neurological Institute (MNI) space (see Figure S2 in the online supplement). Depression severity was assessed with the Quick Inventory of Depressive Symptomatology (QIDS-SR) and the 17-item version of the Hamilton Depression Rating Scale (HAM-D). Although the HAM-D was the primary outcome measure in the THREE-D study, we focused on the QIDS-SR for the results reported in the main text for reasons described in the Methods section of the online supplement. However, we also report effects with the HAM-D in the online supplement. Details on inclusion criteria and rTMS targeting procedures are provided in the Methods section of the online supplement, and the sample characteristics of the 295 THREE-D subjects used in these retrospective analyses are presented in Table 1.

MRI data processing.

All subjects had received a structural brain scan and two 10-minute single-echo resting-state fMRI scans (with the instruction to remain still, stay awake with their eyes closed, and avoid thinking about anything in particular) before and after treatment. MRI data acquisition parameters and scanner type (3T, GE HDx system, 8-channel coil, TR=2 seconds, 32 axial slices, thickness=5 mm) were identical for THREE-D subjects from both sites (Centre for Addiction and Mental Health [CAMH] and University Health Network [UHN]) included in our analyses (19). Acquisition parameters, structural and functional image processing, and denoising procedures are detailed in the Methods section of the online supplement. Unless otherwise specified, r values throughout the Results section represent Pearson correlation coefficients.

Investigating the Impact of Signal Properties on sgACC-StimFC

Effects of common signal-quality metrics (framewise displacement and temporal signal-to-noise ratio [tSNR] [24]) were tested through subsampling across percentiles of each metric. All 590 fMRI scans were independently rated for global signal fluctuations related to burst breathing and deep breaths by three study authors (J.D.P., I.G.E., and C.J.L.) using a 3-point scale indicating the certain presence (2 points), likely presence (1 point), or absence (0 points) of a pattern. The maximum possible score was 12 if all three raters were certain that a particular pattern was present on both scans.

Results

Validation of Imputed sgACC Signal via Weight Maps

One explanation for variable results in studies relating sgACC-StimFC to clinical outcome is that the sgACC signal is degraded and noisy because of its proximity to sinuses in the skull (22). To address this issue, less noisy sgACC signals can be obtained by using a “weight-map” method to impute sgACC signals (2). In brief, this method approximates the sgACC signal by substituting it with the weighted average of all nodes in a scan with time series that are correlated with the sgACC signal in 1,200 subjects from the Human Connectome Project (HCP) (see Figure S1 in the online supplement). We validated this approach in a data set of subjects scanned for up to 15 hours using multiecho fMRI (ME data set), a type of fMRI scan with greatly improved signal qualities in the sgACC (22). In these subjects, individual seed maps of the sgACC were highly reliable across scans (see Figure S3A in the online supplement). This reliability allowed us to construct a “ground truth” sgACC FC map for each subject with their concatenated data sets. The weight-map imputation of the sgACC closely resembled the ground truth seed maps (Figure 1A). In another public data set with 10 30-minute single-echo scans per subject (MSC data set), sgACC seed maps based on the concatenated data (5 hours) were very noisy and were unreliable from scan to scan; in contrast, weight-map approaches with these data produced reliable results in nearly all subjects, with obvious resemblance to the multiecho ground truth (see Figure 1B for representative subjects; see Figure S4 in the online supplement for the full sample). Between-subject variability was preserved with the weight-map method. The effects of sex and breathing patterns on the reliability and accuracy of the weight-map method were also assessed (see the online supplement). Collectively, these results lend strong support to the use of weight maps to impute sgACC signals and to construct reliable sgACC-StimFC maps.

DLPFC Stimulation Modeling Detects Correlation With Treatment Outcomes

Another explanation for variability in relating clinical outcomes to sgACC-StimFC is uncertainty in modeling the stimulated tissue from which the DLPFC signal is composed. We compared two methods of estimating the stimulated DLPFC area in 295 subjects from the THREE-D trial (19) (Table 1). The first method is a heuristic that was often used in earlier studies (2–6, 8, 9) and amounts to a distance-weighted hemisphere (“weighted cone”). The second method involves a biophysical E-field model based on each individual subject’s anatomy and the stimulation parameters. Our primary analysis used the percent improvement in depressive symptoms as measured by the QIDS-SR after completion of the final TMS treatment session as defined in the THREE-D trial. Although the HAM-D was the primary outcome measure in the THREE-D study, our analyses focused on the QIDS-SR (for a rationale, see the “Clinical Outcome Measures” section in the online supplement), and supplemental analyses confirmed similar findings with the HAM-D, as described below.

There was no significant association between sgACC-StimFC and clinical improvement when a generic 12-mm weighted cone was used to approximate the stimulated portion of the DLPFC (Figure 2; see Figure S5 in the online supplement for a range of cone radius thresholds). In contrast, when using individual-specific E-field predictions to define the portion of the DLPFC that was most strongly stimulated (99th percentile of E-field), individual differences in treatment response were correlated with sgACC-StimFC (r=−0.16, p=0.006) (Figure 2). This was true across a range of E-field thresholds (see Figure S5 in the online supplement), in combined (Figure 2) and pretreatment fMRI scans, and after the inclusion of relevant covariates (including treatment condition and sex) with either the QIDS-SR or the HAM-D results as dependent variables (see Results section in the online supplement). No significant results were obtained when using a simple sgACC seed as opposed to the signal imputed from the weight-map approach (see Figure S6A in the online supplement) or when the global cortical signal—removed in previously published positive reports (3, 5, 6, 8)—was not removed (see Figure S6B in the online supplement). Likewise, no association was found when the spatial distance to the most negative sgACC FC site (“optimal site”) was used as a predictor instead of FC strength, which is an approach proposed in recent studies (6, 9) (see Figure S6C in the online supplement).

In summary, a modest association between sgACC-StimFC and treatment outcome was found only when sgACC FC and the stimulated DLPFC area were carefully modeled. For the above analyses, data across treatment conditions and both of a subject’s fMRI scans were collapsed to increase the power to detect a true effect through sample size. This was supported by the absence of rTMS treatment modality effects on sgACC-StimFC (minimal p_FDR=0.29), which was consistent with previous reports (6). However, when analyzed separately, effect sizes were similar across scans (pretreatment scans: r=−0.14; posttreatment scans: r=−0.14), across treatment conditions (iTBS group: r=−0.19; 10 Hz rTMS group: r=−0.13), and across study sites (CAMH sample: r=−0.18; UHN sample: r=−0.13).

To better understand whether the functional connections of the stimulation site to areas other than the sgACC may provide additional treatment prediction, we extended our analysis to the entire brain, using a previously validated multimodal brain parcellation (25). FC of each parcel with the stimulation site was derived by the same weight-map method detailed above for the sgACC analysis. The neuroanatomical distribution and network affiliation of all parcels and the correlation between FC and clinical improvement are depicted in Figure S7 in the online supplement. Seventy parcels showed a nominally significant correlation; the majority of negative correlations mapped onto the default mode network, and the majority of positive correlations mapped onto the cingulo-opercular network and dorsal attention network. Overall, correlations with the left and right sgACC parcels were among the strongest correlations (r=−0.18, rank 4/360, and r=−0.16, rank 18/360, respectively). Taken together, these results indicate that the association of the sgACC with clinical improvement is not a focal property of this brain region but is part of a larger FC signature that involves the default mode network and its anticorrelated networks.

Small Samples Lead to Overestimation and Underestimation of the Whole-Sample Effect

As noted above, previous smaller studies have reported considerably stronger associations between sgACC-StimFC and treatment outcome than the associations reported here. To understand whether sample size limitations might partially explain these conflicting results, we evaluated the impact of sample size on effect size estimates with the THREE-D data set by analyzing the correlation between treatment outcome and sgACC-StimFC in 10,000 random draws for sample sizes ranging from 15 to 280 patients. This analysis revealed substantially greater effect size variance in smaller than larger subsamples, with a wide range of effect sizes at a sample size of 25, which is the median sample size in previous studies (Figure 3A). Even with 50 patients, effects of r=−0.50 or larger were commonly observed. These results suggest that sample sizes alone can account for apparently conflicting results.

sgACC-StimFC Predicts Treatment Outcomes in Patients With Large Global Signal Fluctuations

Another explanation for variability across studies is data quality (24). We tested this by first sorting subjects by the tSNR (mean intensity over the variance in a scan) of their scans. We then computed the effect of sgACC-StimFC on clinical improvement in overlapping subsamples of 50 subjects, going from subjects with the lowest to the highest tSNR values. Unexpectedly, the association between sgACC-StimFC and treatment outcome was strongest in subjects with the lowest tSNR values (Figure 3B). No such effect was found when subjects were sorted by framewise displacement values (see Figure S8 in the online supplement). These results indicate that the effect is accentuated in seemingly lower-quality data. The remainder of this section details our efforts to understand what this subpopulation might represent.

We first disentangled the signal characteristics underlying tSNR in our data and found it to be mainly determined by differences in signal variance (Figure 4A). One major source of variance in fMRI signals stems from common irregular breathing patterns involving deep breaths and burst breathing, which were recently linked to specific blood-oxygen-level-dependent (BOLD) signals (20). While deep breaths cause isolated bands of global BOLD signal fluctuations, burst breathing—characterized by cycles of spindle-like increases and decreases in respiratory depth—produces repeated, serial bands of time-lagged brain-wide BOLD signals. These signals are easily quantifiable when the fMRI time series is visualized as a two-dimensional gray-scaled plot, with cortical nodes as rows and time as columns (“carpet-plot” format). For details on scoring, and a detailed explanation of the carpet-plot format, see Figure S9 in the online supplement and (20). Using the carpet-plot format, three blinded raters (J.D.P., C.J.L., and I.G.E.) scored the BOLD signal time series of all 295 subjects (590 fMRI scans) for global signal patterns indicative of burst breathing and deep breaths (Figure 4B). Ratings were highly similar across raters (Cohen’s kappa values ranged from 0.59 to 0.75) (Figure 4C, left panel) and replicated previously observed sex biases (20) (Figure 4C, right panel). Together, these results validate our approach to identifying global signal signatures that have been previously linked to two forms of irregular breathing.

Comparisons between subjects with strong evidence for burst breathing (N=47) relative to the remainder of the sample (N=248) revealed several interesting observations. First, global signal patterns indicative of burst breathing were associated with lower tSNR, as hypothesized, given the increased BOLD signal variance induced by bursting (see Figure S10 in the online supplement). Second, among the subjects displaying global signal patterns related to burst breathing, sgACC-StimFC was strongly correlated with clinical improvement (r=−0.49, p=0.0004) (Figure 4D, left panel) (see Figure S11 in the online supplement for a range of cutoff values to assess burst certainty). The probability of drawing a subsample of that size with an effect this strong by chance is low (p=0.002 from 10,000 samples selected at random without replacement) (Figure 4D, right panel). Moreover, randomly subsampling subjects across percentiles of bursting scores showed the effect to increase as a function of a sample’s bursting score (see Figure S12 in the online supplement). Third, when the subsample of 47 subjects displaying global signal fluctuations reflecting burst breathing was removed, the association between sgACC-StimFC and clinical improvement was clearly weakened in the remaining sample (N=248, r=−0.09, p=0.32) (Figure 4E, left panel). These effects were specific to subjects with evidence of bursting and not a common property of subjects with low tSNR scans: a control group of 47 subjects with equally low tSNR (105.56 vs. 117.48 in those showing burst breathing) but no evidence of burst breathing did not show an accentuated effect (N=47, r=−0.12, p=0.42) (Figure 4E, right panel). Lastly, the dependency of the effect on tSNR was removed with the absence of the 47 subjects showing evidence of burst breathing (Figure 4F), indicating that this phenomenon was driven by the presence of burst breathing in the sample rather than being a consequence of nonspecific sources of high signal variance. Taken together, these results indicate that the link between sgACC-StimFC and clinical improvement in our sample is explained by a specific subset of subjects with global signal fluctuations that have previously been linked to burst breathing.

Discussion

Here, in a large rTMS data set, we found a relationship between sgACC-StimFC and treatment outcome, as previously reported. The effect we observed was substantially smaller than that described in several previous reports and accounted for only 3% of treatment response variance, hence calling into question the utility of this marker in its current form. Furthermore, this relationship was obtained only when weight maps were used to impute otherwise noisy sgACC FC, when the stimulated part of the DLPFC was estimated with E-field modeling, and when the sample included certain subjects, namely, those with low tSNR due to high signal fluctuations that have been specifically linked to irregular breathing. The interpretation of these findings is multifaceted. On the one hand, the association between sgACC-StimFC and treatment outcomes was modest in this sample and was contingent on certain choices. On the other hand, the effect was strengthened by more accurate modeling approaches, and it was robust to the partitioning of the data, robust across measurement scales, and robust to the inclusion of covariates. Furthermore, it was concentrated in a particular set of subjects who display identifiable characteristics. Our interpretation is at once tempered and hopeful: striking, strong linkages of FC to outcome in a large sample were not found, but the concentration of these linkages in a subpopulation suggests potential avenues toward understanding who is responsive to treatment and why.

Unexpectedly, the predictive value of sgACC-StimFC was strongest in patients with patterns of global signal variance associated with burst breathing. One possible explanation derives from the recent observation that burst breathing coincides with high-amplitude BOLD events that are time-lagged between networks (20). The spatiotemporal structure of these time-lagged events is clearly revealed after removing CO₂-related global signal fluctuations and is such that it accentuates negative FC between the default mode network (of which the sgACC is commonly a part) and networks represented in the DLPFC (salience and frontoparietal networks). Consistent with this notion, subjects with prominent burst breathing had more negative sgACC-StimFC (see Figure S10B in the online supplement), and omission of global signal regression removed the observed effect (see Figure S6B in the online supplement). Hence, the anticorrelation in FC between the stimulation site and the sgACC might be an indirect index of what intrinsic networks have been engaged by TMS, which might be the true mediator of the effect. This would be aligned with a rich literature linking individually variable functional network topology to various behavioral domains (21, 26, 27). Future work might elucidate the utility of these more comprehensive measures of functional brain organization that are only obtainable from high-quality data, such as multiecho fMRI scans. Alternatively, a propensity for burst-like respiration may differentiate a subsample of depressed patients for whom sgACC FC strongly mediates rTMS outcomes. We did not find strong evidence in the MSC and ME data for burst breathing having an impact on the validity of the weight-map method, besides a minimally better run-to-run reliability in male versus female MSC subjects (which could be indirect evidence, as burst breathing is more prevalent in males; see Results section in the online supplement). If reliability is increased in subjects with burst breathing, this could increase the power to detect a true effect in these subjects. All of these factors may contribute to the effect. It is also possible that burst breathing may be influenced by somnolence or fluctuating levels of arousal during the scan. However, at present, there is no evidence linking burst breathing to low vigilance or somnolence, and several observations in the HCP data speak against it (20), rendering vigilance an unlikely mediator of these effects.

With regard to data quality, it should be noted that the mean tSNR values were high for BOLD fMRI data, which was likely due to the somewhat larger-than-average voxel size. In terms of acquisition parameters, the THREE-D data compare favorably to data from previous studies (5, 6, 9, 16–18), with the shortest TR (2 seconds) and second-longest acquisition time (total of 20 minutes) among those studies. These factors render data quality an unlikely candidate to explain effect size discrepancies between our findings and those of previous studies.

We discuss other sources of heterogeneity that potentially contribute to the discrepant effect sizes between those of our study and previous studies in the Discussion section in the online supplement. One possible source of heterogeneity is the study population: major depressive disorder is not a unitary disease entity, and sgACC-StimFC might be more relevant for predicting rTMS outcome in particular subpopulations with major depressive disorder. The THREE-D sample comprised patients with major depressive disorder and a relatively high rate of comorbid anxiety and polypharmacy (19), as seen in naturalistic samples. However, previous reports of large effect sizes also included populations with similarly high rates of comorbid anxiety and polypharmacy (5, 6), indicating that these variables alone are unlikely to explain the small effect we observed. The effect size we report (r=−0.16) represents a “best-case scenario” in our sample, where data across scans and treatment conditions were concatenated to increase power, and results with the QIDS-SR, which yielded stronger effects than the HAM-D, were used. Consequently, the small effect size might itself be inflated and likely represents the upper bound of a true effect obtainable in this data set.

Several limitations should be noted. First, like all previous studies, this was a retrospective analysis as opposed to a prospective analysis with a preformulated analysis pipeline. Second, subjects in the THREE-D study received either 10 Hz rTMS or iTBS, while previous studies of sgACC FC investigated the effects with 10 Hz rTMS only. However, we found a similar effect in both conditions, which suggests that our conclusions apply to both 10 Hz rTMS and iTBS treatment modalities. Third, like other studies of sgACC FC and rTMS, the THREE-D study did not include a sham control condition. Both the sgACC and DLPFC have been implicated in placebo responses (28, 29), and our observation of a modest correlation between symptom improvement and sgACC-StimFC might be attributable to treatment responses, placebo responses, or both. The DLPFC, in particular, has shown consistent involvement in placebo antidepressant responses across metabolic, perfusion-based, and fMRI studies (30). On the other hand, we also found that accurate E-field modeling strengthened the association between sgACC-StimFC and outcome, which may point to a treatment effect. E-field modeling of the stimulation site decreases variance by sampling FC features that were stimulated with TMS but increases subject-to-subject variance with respect to which FC features are being averaged within the DLPFC (i.e., the chosen features depend on each subject’s unique E-field distribution and shape). Since the neuroanatomical substrate of a placebo response is unlikely to move from subject to subject in line with their E-field, we interpret this as supporting evidence for a TMS treatment effect and not just a placebo response. Fourth, neither our study nor earlier ones included pupillometric or electrophysiological monitoring of vigilance, which can have marked confounding effects on the global signal and FC. A further limitation was the lack of respiratory belt traces in the THREE-D data. However, bursting signal patterns closely resembled previously reported patterns in the HCP data (Figure 4B), and we exactly replicated a burst breathing–specific sex bias previously reported in the HCP data (20), which is supportive evidence for the correct identification of irregular breathing patterns in the THREE-D data set (Figure 4C). Lastly, the validity of our results—and those of previous studies—rests on the accurate identification of the stimulated cortex, which is an active area of research (23). We used a biophysical E-field model that is, to our knowledge, the primary method for estimating the spatial distribution of TMS effects. However, the extent to which this model accurately predicts neuronal populations modulated by rTMS is still partially unknown, and future models might yield better approximations. Additional discussion of alternative explanations of our findings associated with burst breathing, our effect size estimate, and analytical choices is provided in the Discussion section in the online supplement.

Together, these results suggest the need to reevaluate the utility of current sgACC FC-based targeting approaches for rTMS and elucidate relevant sources of variability that might have led to mixed results in previous studies. At the same time, these results provide strong evidence for the existence of a true—albeit weak—association. This highlights the need to further explore fMRI predictors for rTMS targeting, possibly involving additional network measures and fMRI data types. Methodological improvements, including multiecho fMRI and dense-sampling approaches, have already shown promise for enabling more comprehensive and personalized characterizations of functional brain organization that could be used to make better treatment predictions in the future. Prospective clinical trials designed to evaluate the clinical utility of fMRI-guided targeting strategies in naturalistic, real-world settings will be critical.