Anish Mitra,Marcus E. Raichle,Andrew D. Geoly,Ian H. Kratter,Nolan R. Williams
Targeted neurostimulation reverses a spatiotemporal biomarker of treatment-resistant depression
depression,brain stimulation,functional MRI
Significance Major depressive disorder (MDD) is one of the most significant causes of morbidity in the world. However, despite decades of functional neuroimaging research, the biological signatures of MDD remain poorly understood. Here, we characterize the propagation of spontaneous brain activity using fMRI to investigate whether the flow of neural activity differs between patients with MDD and typical controls. We find a group of patients with MDD whose symptoms are linked to aberrant propagation of brain activity within the salience network. Additionally, TMS can restore typical propagation patterns in the salience network and relieve symptoms. In sum, our results offer a pathway for identifying and treating mechanisms of neuropsychiatric diseases.
Targeted neurostimulation reverses a spatiotemporal biomarker of treatment-resistant depression

Major depressive disorder (MDD) is one of the most significant causes of morbidity in the world. However, despite decades of functional neuroimaging research, the biological signatures of MDD remain poorly understood. Here, we characterize the propagation of spontaneous brain activity using fMRI to investigate whether the flow of neural activity differs between patients with MDD and typical controls. We find a group of patients with MDD whose symptoms are linked to aberrant propagation of brain activity within the salience network. Additionally, TMS can restore typical propagation patterns in the salience network and relieve symptoms. In sum, our results offer a pathway for identifying and treating mechanisms of neuropsychiatric diseases.

Major depressive disorder (MDD) is the leading cause of disability worldwide (1). Human neuroimaging has revealed that the neural mechanisms of MDD are not restricted to any single brain region, but, instead, likely involve aberrant function across brain-wide networks (2–14). Several studies have specifically postulated that imbalanced signaling between the anterior cingulate cortex (ACC) and other brain areas involved in emotional processing may underlie MDD (2, 7, 15–19). Indeed, repetitive transcranial magnetic stimulation (rTMS) protocols have been theorized to modulate signaling between the ACC and other brain regions by stimulating the left dorsolateral prefrontal cortex (DLPFC) (2, 7, 16, 20–27). Yet, progress in developing treatments for MDD has been hindered by our inability to reliably detect mechanistic biomarkers to inform diagnosis, treatment, and development of new therapeutics. The search for neurobiological biomarkers of MDD has been slowed by both conceptual and experimental barriers. First, although several resting-state-functional MRI (rs-fMRI) studies of human brain function have reported functional differences between patients with MDD and healthy controls (HCs), these differences are rarely subject to targeted network manipulations. Consequently, it has been difficult to assign causal relationships between mechanistic biomarkers and MDD. Second, although MDD is believed to arise from disordered communication in the brain, the vast majority of rs-fMRI studies examine either activity levels within individual brain regions or interactions between brain regions using zero-lag correlation (conventionally called “functional connectivity”) (2, 4, 6, 8, 10, 22, 23, 28). A limitation of conventional functional connectivity is that there is no assignment of signaling directionality. Thus, a circumstance in which brain region A signals to brain region B can appear identical to a circumstance in which the direction of communication is reversed, even though the neurobiological implications of signaling direction may be profound. Recent work has found that rs-fMRI contains a robust, highly reproducible temporal structure characterized by stereotyped directional signaling patterns across the human brain (29–34). Directed patterns in rs-fMRI correspond with infraslow (<0.01 Hz) human electrocorticography, mouse optical calcium imaging, and mouse electrophysiology (35, 36). Additionally, there is evidence that the directional flow of rs-fMRI activity is exquisitely sensitive to brain state, more so than the zero-lag correlations that define conventional functional connectivity (36, 37). Taken together, these developments raise the intriguing possibility that directional signaling patterns in rs-fMRI may be a sensitive marker of aberrant network activity in neuropsychiatric disorders such as MDD. Here, we investigate the relationship between brain-wide directed signaling patterns and MDD by leveraging rs-fMRI data acquired during clinical trials for our FDA-cleared neuromodulation approach. In brief, our MDD therapeutic approach, termed Stanford neuromodulation therapy (SNT), previously published in two open-label studies as well as a randomized controlled trial (16, 20, 21), involves daily administration of 10 neurostimulation sessions (applied through an rTMS device) each day for 5 d using a spaced intermittent theta burst pulse design applied to individualized left DLPFC targets (20, 21). While the SNT depression protocol has proven effective in the acute clinical treatment of treatment-resistant MDD (TRD) (16, 20, 21), the biological basis of the treatment effect remains unknown. In the present study, we analyze directed signaling patterns in rs-fMRI data acquired immediately before and after SNT treatment in two clinical trials: First, a double-blind randomized clinical trial in which patients with treatment resistant depression (TRD) received either active or sham SNT, and second, an open-label clinical trial in which patients with TRD received active SNT. In addition, we also analyze brain-wide directed signaling patterns in typical control subjects without MDD (Table 1). By combining these datasets, we test: 1) whether SNT modulates directed signaling in rs-fMRI (beyond placebo effects in sham SNT), 2) whether SNT-induced modulation of directed signaling in rs-fMRI relates to symptomatic improvement of MDD, 3) whether directed signaling patterns in TRD differ from typical controls, 4) whether SNT restores a typical flow of directed activity in rs-fMRI, and 5) whether baseline rs-fMRI directed signaling patterns predict SNT treatment response. To preview our results, we find that SNT modulates the directional flow of rs-fMRI activity between the dorsal ACC and the salience network, which has previously been implicated in emotional processing (8, 41, 42). Variance in the magnitude of the temporal effect corresponds with treatment efficacy. Furthermore, we find that SNT corrects a baseline aberrancy in directed rs-fMRI signals and that baseline rs-fMRI temporal structure predicts both individual depression severity and treatment response.

Our method for computing the directed activity in rs-fMRI data has been described in detail in prior publications and is briefly summarized in Fig. 1 (see Methods and refs. (31, 33, 35) for further details). Fig. 1A shows two exemplar rs-fMRI time series from two brain areas. On visual inspection, it is clear that the blue time series generally leads the red time series. Yet, the temporal difference between the two time series is small, less than the two second sampling density of this rs-fMRI example. We compute the temporal delay between these time series by computing the lagged correlation function (Fig. 1B). The peak of the lagged correlation function represents the temporal delay corresponding to the average direction and temporal offset between the two time series (31). Although the empirical peak is at zero-lag in Fig. 1B (yellow dots, black line), it is apparent that the lagged correlation function is asymmetric about the origin, indicating the existence of a non-zero temporal delay. We recover this delay using parabolic interpolation (green line in Fig. 1B), which indicates a nonzero peak to the lagged correlation function (orange dot in Fig. 1B) occurring at a temporal offset of approximately negative 0.6 s, indicating that the blue time series leads the red time series. Note that this delay does not mean that the blue time series always leads the red time series by 0.6 s (there are evident counterexamples in Fig. 1A). Instead, the computed delay reflects an average over the whole time series. Our subsequent temporal analyses of rs-fMRI data are computed using parabolic interpolation on lagged correlation curves. Prior work has validated that the interpolated delays calculated in rs-fMRI data reflect similar delays in underlying neural activity (35, 36).

We first test the hypothesis that SNT induces changes in the temporal structure of resting-state-fMRI not found in sham SNT by analyzing rs-fMRI data collected as part of a double-blind, randomized controlled trial of SNT in outpatients with TRD (Dataset 1, Table 1). In the trial, patients were randomized to receive 5 d of either active SNT or a sham control (which utilized the same coil with no outward indication of sham vs. active treatment; see Methods) over their left DLPFC. rs-fMRI data and behavioral metrics were acquired at baseline prior to SNT treatment (“baseline” in Fig. 1) and again posttreatment (“1 wk” in Fig. 1) within 3 d of treatment completion. Of note, the present analysis tests whether posttreatment directed signaling patterns differ from baseline; we did not acquire fMRI data concurrent with rTMS application. To analyze directed signaling patterns, we compute lags for all pairs of voxels in gray matter in rs-fMRI data (Fig. 1 A and B). These results are assembled into time-delay (TD) matrices, which have dimensions voxels × voxels and entries in units of seconds (Fig. 1C). TD matrices represent the lag between all pairs of voxels in gray matter. We then compute the mean over all columns of the TD matrix to yield a lag projection map (Fig. 1C). Lag projection maps topographically represent the mean lag between each voxel and the rest of the brain. As we are interested in treatment-induced changes, we compute the lag projection at baseline and posttreatment. We then compute a lag projection difference map for each subject by simply subtracting the baseline lag projection from the 1-wk lag projection. The lag projection difference map captures subject-wise changes in brain-wide directed rs-fMRI signaling. To identify statistically significant spatial clusters attributable to active SNT, we next compute threshold–extent criteria by extensive permutation resampling (Methods) comparing active SNT lag projection difference maps vs. sham SNT lag projection difference maps. The outcome of this comparison identifies brain areas with significant temporal shifts as a result of active SNT treatment that are not found in the sham-SNT group. The result, shown in Fig. 1D, shows that there are two regions with statistically significant temporal shifts attributable to active SNT: the DLPFC underlying the SNT stimulation site (SI Appendix, Fig. S1) and the bilateral ACC. The left DLPFC (cool cues) shifts earlier with respect to the rest of the brain as a function of the SNT protocol, whereas the ACC (warm hues) shifts later with respect to the rest of the brain. Having found two voxel clusters exhibiting statistically significant changes in rs-fMRI temporal structure in the whole-brain active-SNT vs. sham-SNT contrast, we next compute seed-based lag maps using the clusters shown in Fig. 1D as seeds (Fig. 1C, Methods) to examine the specific network relationships driving the brain-wide effects. Seed-based lag maps represent temporal lags between each voxel and the average time course computed over the seed-region of interest. Thus, whereas lag projections reveal the average temporal position of a voxel with respect to the rest of the brain, seed-based lag maps depict pair-wise temporal relationships between the seed-region and every other voxel. Cool colors depict regions that are earlier than the seed and warm colors depict regions that are later than the seed. Fig. 2 displays the group-average seed-based lag map for the left DLPFC and ACC seeds at both baseline and posttreatment in active SNT. Statistically significant differences across baseline and the posttreatment seed-based lag maps, computed by permutation resampling, are shown in the far-right column of Fig. 2. Note that the statistical comparison of seed-based lag maps in Fig. 2 A and B derived from Fig. 1D represents a post-hoc analysis to understand which specific pair-wise effects most reliably contribute to the observed differences in lag projections. As evident in Fig. 2A, there is a statistically significant temporal shift between the left DLPFC and ACC in active SNT, such that the ACC becomes later than the left DLPFC. Thus, the increased “earliness” of the left DLPFC in the brain-wide analysis of Fig. 1D is driven most strongly by its relative shift with respect to ACC. Fig. 2B depicts that there is a post-hoc statistically significant shift between the ACC and several regions, including the anterior insula, lateral prefrontal cortex, and temporoparietal junction. Each of these regions shifts earlier with respect to the ACC as a function of SNT, thus driving the increased brain-wide “lateness” of ACC depicted in Fig. 1C. Strikingly, several areas, including the anterior insula and lateral prefrontal cortex, reverse their temporal ordering with respect to the ACC. At baseline, in patients with TRD, anterior insula and lateral prefrontal cortex are late with respect to the ACC (Fig. 2B, “Baseline”). However, after SNT treatment, these structures shift to becoming earlier than the ACC (Fig. 2B, “1 wk”). To test the reproducibility of SNT-induced temporal shifts, we compute statistically significant lag map effects for the left DLPFC and bilateral ACC clusters derived in Fig. 1D in a separate, open-label SNT study (Dataset 2, Table 1). The open-label study was identical to the randomized clinical trial except all subjects in the open-label study received active SNT and were aware of this fact. Despite this difference, our analysis of left DLPFC-seeded lag maps shows that the ACC exhibits a statistically significant shift later than that of the left DLPFC (Fig. 3A) as a function of SNT, mirroring the results from our randomized clinical trial. Similarly, our analysis of ACC-seeded lag maps (Fig. 3B) shows that anterior insula, lateral prefrontal cortex, and temporoparietal junction shift earlier with respect to the ACC as a function of our SNT protocol. Finally, to test whether the temporal shifts we observed are accompanied by changes in correlation, we compute zero-lag correlations for both the left DLPFC and ACC across baseline and posttreatment rs-fMRI scans. The results, shown in SI Appendix, Figs. S2 and S3, reveal no significant changes in ACC- or left DLPFC-seeded correlation structure as a function of SNT. Thus, active SNT induces focal shifts in rs-fMRI temporal structure in the left DLPFC and bilateral ACC. The left DLPFC moves earlier in our brain-wide analysis, with a statistically significant shift with respect to the ACC. The bilateral ACC moves later in our brain-wide analysis, with statistically significant shifts with respect to the anterior insula, lateral prefrontal cortex, and temporoparietal junction. We find highly similar results in a second, independent open-label trial of SNT.

Having found focal SNT-induced temporal shifts in the left DLPFC and ACC, we next ask whether variability in temporal shifts in these brain regions is associated with treatment response. We explore this question by computing seed-based left DLPFC and ACC temporal shifts in each subject across both the randomized clinical trial (Dataset 1) and the open-label study (Dataset 2), restricted to a conjunction of the statistically significant spatial clusters found in Figs. 2 and 3 A and B. Critically, the same seeds are used in all subjects across both datasets (Seed-Mask Construction in Methods). Thus, for the left DLPFC seed, we compute in each individual the temporal shift between the left DLPFC and a conjunction of the anterior cingulate clusters shown in Figs. 2 and 3A. Similarly, for the ACC seed, we compute in each individual the temporal shift between the ACC and a conjunction of the anterior insula, lateral prefrontal cortex, and temporoparietal junction clusters shown in Figs. 2 and 3B. Fig. 3C shows the relationship between the left DLPFC temporal shift and improvement in depression symptoms (quantified by percent improvement in MADRS score) across individuals. As dictated by our prior results, the left DLPFC shifts early in the vast majority of subjects over 1 wk of treatment. However, there is no statistically significant relationship between the extent of the left DLPFC shift toward earliness and improvement in depression symptoms. In contrast, Fig. 3D shows a robust, statistically significant relationship between a late shift in the ACC and improvement in depression symptoms (Spearman’s , corrected). The relationship between the ACC temporal shift and improvement in MADRS scores is logarithmic (solid line in Fig. 3D), suggesting that there are saturation effects to benefits of the late shift in ACC. We also computed the relationship between the ACC temporal shift and improvement in MADRS scores in subjects given sham treatment and observed a weaker nonsignificant effect (SI Appendix, Fig. S4). Thus, whereas SNT induces both an early shift by the left DLPFC and a late shift by the ACC, variability in treatment response relates specifically to the late shift by the ACC with respect to the anterior insula, lateral prefrontal cortex, and temporoparietal junction.

Given the relationship between the SNT-induced temporal shift in ACC signaling and improvement in depression symptoms, we next test the hypothesis that aberrant ACC signaling is a biomarker of MDD that is reversed by the SNT protocol. We approach this question by computing seed-based lag maps of the ACC and comparing across patients with major depression and HCs (Seed-Map Construction in Methods). For our TRD cohort, we include the baseline rs-fMRI scans across both the randomized clinical trial and the open-label study (Dataset 1 and Dataset 2; N = 33). For comparison, we use two HC cohorts. rs-fMRI data in the first group of HCs (HC1; Dataset 3 in Table 1) were collected on the same scanner as our TRD cohort using the same scan parameters; however, the sample size is relatively modest (N = 16). rs-fMRI data in the second group of HCs (HC2; Dataset 4 in Table 1) were collected on a different fMRI scanner using different scan parameters (Methods) but in a larger group of subjects (N = 85). Given the importance of reproducibility in fMRI, we intentionally chose to test the consistency of our findings in a large set of rs-fMRI data acquired using different parameters than data acquired in-house. The ACC-seeded lag maps for each group are shown in Fig. 4 A–C. Fig. 4 D and E depict statistically significant differences between the TRD cohort and the two HC cohorts (HC1 and HC2), respectively. The most prominent effect in Fig. 4 D and E is found in the anterior insula bilaterally, which is late with respect to the ACC in our TRD cohort, but early with respect to the ACC in both of our HC cohorts. Further effects include lateral prefrontal cortex bilaterally, which is later with respect to the ACC in TRD than that in our HC cohorts. Therefore, our results show that directed rs-fMRI signaling through the ACC distinguishes patients from HCs. Specifically, the anterior insula is late with respect to the ACC in MDD but early with respect to the ACC in HCs, whereas the lateral prefrontal cortex is later with respect to the ACC in MDD than that in HCs. Given that the SNT protocol shifts these regions earlier with respect to the ACC (Fig. 2B), our results imply that SNT shifts the temporal structure of ACC-seeded rs-fMRI activity from an aberrant state in MDD toward a state that resembles HCs. However, comparing Fig. 2B (“Post-Treat”) with Fig. 4 B and C, it is also apparent that SNT also induces a nonphysiologic change in the bilateral DLPFC, such that the DLPFC becomes earlier than the ACC post-SNT but is later than the ACC in HCs.

In light of the observed clinical variability in response to the SNT protocol, we next test whether the pretreatment ACC-based signaling in rs-fMRI activity predicts treatment response. In order to visualize potential differences, we computed ACC-seeded lag maps (Seed-Map Construction in Methods) in subjects with high response to the SNT protocol (>90% reduction in MADRS score, N = 6) against subjects who were nonresponders to the SNT protocol (<50% reduction in MADRS score, N = 6). The results, illustrated in Fig. 5, reveal that the “high-response” group (Fig. 5A) exhibit an ACC-seeded lag map similar to the TRD group average shown in Fig. 4A. In particular, note that the anterior insula is late with respect to the ACC, and that the lateral prefrontal cortex is quite late with respect to the ACC (~0.5 s). In contrast, the “nonresponse” group (Fig. 5B) resembles the HC group averages in Fig. 4 B and C. Specifically, the anterior insula is early with respect to the ACC, and prefrontal cortices are far less late with respect to the ACC (~0.2 s). Statistical differences in ACC-seeded lag maps between the “high-response” and “nonresponse” groups are shown in Fig. 4C, which again illustrates that the anterior insula and lateral prefrontal cortex are earlier with respect to the ACC in the “nonresponse” group as compared to the “high-response” group. Although the differences between the “high-response” and “low-response” groups visualized in Fig. 5 A–C suggest a general relationship, we must explicitly test this hypothesis across all subjects. Thus, we next examine whether baseline ACC temporal structure is linked to treatment response in all 23 patients who received the SNT protocol. The result, shown in Fig. 5D, reveals that there is indeed a statistically significant relationship between ACC temporal structure at pretreatment baseline and subsequent SNT treatment response (lags between the ACC and network partners computed as in Fig. 3D, see Seed-Based map Construction in Methods). Thus, patients with low response to the SNT protocol exhibit baseline ACC-seeded lag maps that resemble HCs, whereas patients with high response exhibit highly aberrant ACC-seeded lag maps.

In analyzing the relationship between baseline ACC temporal structure in the “high-response” and “nonresponse” groups, we found that the two groups were not matched in terms of depression severity as measured by the MADRS score. Instead, the “high-response” group exhibited higher baseline MADRS scores (more severe depression) as compared to the “nonresponse” group. We thus test whether ACC temporal structure is also related to depression severity across all 33 patients in the present study (active, sham, and open label). The result, shown in Fig. 5E, illustrates that there is a statistically significant relationship between the baseline ACC temporal structure and MADRS score in our TRD patients.

The present study investigates the relationship between directed signaling patterns in human rs-fMRI and MDD both in assessing differences between patients and controls and the effects of SNT, a recently developed rTMS protocol. In a double-blind randomized clinical trial, we found that active SNT (but not sham SNT) caused specific directed-signaling shifts in the left DLPFC and bilateral ACC (Figs. 1D and 2). We found highly similar effects in a separate open-label trial (Fig. 3 A and B). Whereas the temporal shift of the left DLPFC is not related to treatment efficacy, the temporal shift in the ACC with respect to brain regions including the anterior insula, lateral prefrontal cortex, and temporoparietal junction is correlated with symptomatic improvement (Fig. 3 C and D). Comparing between patients and HCs, we further found that the ACC is abnormally early with respect to the anterior insula and DLPFC in MDD (Fig. 4), indicating that SNT reverses, in part, aberrant directional signaling patterns. Finally, we examined associations between baseline depression severity, treatment response, and directed signaling to find that aberrant ACC-earliness is correlated with both greater baseline MADRS scores and greater treatment response to SNT (Fig. 5). Taken together, our findings suggest spatiotemporal biomarker of MDD on the basis of aberrant directed signaling between the ACC, anterior insula, lateral prefrontal cortex, and temporoparietal junction. The combination of these brain areas is a well-known brain network that has been called the salience network (42), the cingulo-opercular network (43), the ventral attention network (44), or the midcingulo-insular network (45), among other monikers. For simplicity, we shall refer to this constellation of brain areas as the salience network. Our findings are in conceptual accordance with a large body of task-based fMRI studies that have linked activity within the salience network to emotional processing and evaluations of internal state (10, 41, 46–50). Moreover, at least one prior study has found that fMRI activity in the anterior insula leads the ACC during a salience task, consistent with our resting-state finding that in HCs, the anterior insula leads the ACC (51). Conventional (zero-lag) functional connectivity differences within the salience network have also been reported in MDD and in response to rTMS applied to the left DLPFC, although the direction of functional connectivity change (increase or decrease) is variable (8, 17, 24, 52). Our analysis of conventional functional connectivity of the left DLPFC and ACC found no statistical differences between patients and controls nor between pretreatment and posttreatment scans (SI Appendix, Figs. S2 and S3). This discrepancy may be attributable to methodological differences in quantifying changes in functional connectivity. However, our findings are consistent with recent work showing that there is a weak relationship between conventional connectivity and behavioral measures (53) and add to growing evidence that directional rs-fMRI signaling may be a more potent marker of brain-state than zero-lag correlations (36, 37). Notably, our results highlight the key functional role of a major hub in the salience network, the dorsal ACC. Prior reports have variably linked the pathophysiology underlying MDD as well as the mechanism of treatment of rTMS applied to the left DLPFC to both dorsal ACC (4, 9, 23, 28, 54–56) and subgenual ACC (5, 7, 23, 28, 57). Indeed, the left DLPFC target in SNT is derived by identifying the portion of the left DLPFC is most anti-correlated with subgenual ACC. However, both subgenual and dorsal ACC are anticorrelated with the left DLPFC [SI Appendix, Fig. S2 (22); thus, modulation of the left DLPFC has the potential to alter activity in both regions of the ACC. Future hypothesis-directed studies are needed to specifically test whether directional signaling is altered in subgenual ACC and other brain regions that have been previously implicated in MDD. On the basis of our findings, we propose a model of an MDD biotype in which the pathophysiology is marked by aberrantly early signaling by the dorsal ACC within the salience network. However, several questions remain. First, given work highlighting successful treatment of MDD through alternative brain regions (14, 57, 58), further studies must determine whether the dorsal ACC-salience mechanism found here is a unique biotype of MDD or a part of an overlapping set of neurophysiologic causes of MDD (6). Second, our hypothesis regarding the mechanistic function of SNT derives from a recently developed treatment paradigm involving a 5-d course of spaced intermittent theta-burst stimulation (see Methods, and refs. 20 and 21 for further details). Future work is required to determine the role of administration frequency and pulse design in how rTMS modulates the temporal structure of rs-fMRI. Third, our patient cohort consists of individuals with treatment-resistant depression. Moreover, the biomarker we have identified is associated with more severe depression (higher MADRS scores) in our sample (Fig. 5E). Thus, future studies are needed to understand the generalizability to both treatment-naïve and less severe MDD. Finally, the importance of the temporal relationship between dorsal ACC and anterior insula in our results raises the intriguing possibility that modulation of the anterior insula may be an alternative means for normalizing the temporal flow of activity in the salience network. Recent work suggests that it is possible to noninvasively modulate insular activity (59, 60) (but see also: ref. 61). However, future investigations are required to determine whether modulation of the anterior insula can be utilized for antidepressant effect.

All procedures were conducted in accordance with the ethical standards outlined in the Declaration of Helsinki. The study was approved by the Stanford University Institutional Review Board, and all participants provided written consent before taking part in any study procedures. Four datasets were analyzed in the present report. In dataset 1 (Table 1), we aimed to acquire baseline and posttreatment rs-fMRI data in 29 patients enrolled in our previously published double-blind randomized control trial (20), which contains a detailed description of the trial and patient characteristics. Briefly, each of the patients in dataset 1 had a diagnosis of TRD, where diagnosis of MDD was confirmed with the Mini International Neuropsychiatric Interview (MINI). Participants were required to have not responded to at least one antidepressant medication (minimum trial duration of 6 wk) to be eligible for the study. The participants had no history of a psychotic disorder, substance use disorder, major systemic illness, current mania (assessed using the Young Mania Rating Scale), lesional neurological disorder, or any contraindications to rTMS or MRI (such as implanted metal or history of seizures). Our protocol called for each patient to undergo a baseline MRI scan within 48 h of starting the clinical trial, and then to undergo a second posttreatment MRI scan within 48 h of completing the trial. Out of the 29 patients enrolled in the trial, a total of 24 patients completed both baseline and posttreatment MRI scans. Out of these 24 subjects, 4 were excluded due to excessive head motion in the scanner (see below for criteria; 2 patients were excluded from the active and sham arms, respectively). Thus, there remained 10 subjects each in the “active” and “sham” SNT treatment fMRI analyses. Dataset 2 (Table 1) was collected under identical conditions as Data 1 with the exception that the patients were recruited in an open-label trial, clinical data from which have also been previously reported (21). We enrolled a total of 20 rTMS naïve patients in the trial, 16 of whom completed MRI scans pre- and post-treatment as described previously. Out of these 16 subjects, 3 patients were excluded due to excessive head motion in the scanner (see below for criteria). Dataset 3 (Table 1) was similarly collected under identical conditions to Datasets 1 and 2 with the exception that instead of patients with MDD, we recruited HCs that were demographically matched to our patient samples. HC subjects did not endorse current symptoms of MDD as assessed by MADRS scores and the MINI. Moreover, as in our patient population, our HCs had no history of a psychotic disorder, substance use disorder, major systemic illness, current mania (assessed using the Young Mania Rating Scale), lesional neurological disorder, or any contraindications to rTMS or MRI (such as implanted metal or history of seizures). As opposed to Datasets 1–3, which were all acquired at the same scanner at Stanford with common inclusion/exclusion criteria, Dataset 4 (Table 1) was obtained from the publicly available Harvard-MGH Brain Genomics Superstruct Project (38–40). Subjects in the Superstruct project denied any current/past history of Axis I pathology or neurological disorder, current psychotropic medication usage and/or acute physical illness, or displayed atypical brain anatomy. Moreover, we only included subjects with absolutely minimal motion (fewer than 5% of frames censored due to high motion) and subjects for whom profile of mood state data was available and showed total mood disturbance scores of <0. There were 85 subjects who met each of these criteria.

Participants in Datasets 1 and 2 received a 5-d treatment course under the SNT protocol. Stimulation location was chosen based on an algorithm that selected the area within the DLPFC with the most negative correlation with an anatomically defined subgenual cingulate seed. Further details are contained in our prior reports (20, 21). A MagVenture MagPro X100 (MagVenture A/S, Denmark) system was used to deliver sessions of intermittent thetaburst stimulation: 60 cycles of 10 bursts of 3 pulses at 50 Hz were delivered in 2 s trains (5 Hz) with an 8-s intertrain interval. Stimulation sessions were delivered hourly. Ten sessions were applied per day (18,000 pulses/day) for 5 consecutive days (90,000 pulses in total). Stimulation was delivered at 90% resting motor threshold (rMT). A depth correction (20, 21) was applied to the rMT to adjust for differences in the cortical depth of the individual’s functional target compared to the primary motor cortex in order to consistently achieve 90% rMT in the intended functional target, but stimulation was never delivered above 120% rMT for safety. A Localite Neuronavigation System (Localite GmbH, Sankt Augustin, Germany) was used to position the TMS coil over the individualized functional stimulation target. All participants in this report completed the full 5-d treatment course.

Participants’ depressive symptoms were assessed using the Montgomery–Asberg Depression Rating Scale (MADRS) (62). MADRS assessments were performed before SNT (baseline) and after the 5-d SNT treatment (posttreatment).

All MRI data were acquired either at Stanford (Datasets 1–3, Table 1) or at MGH (Dataset 4) (40). The MGH data were acquired on a 3T Siemens Tim Trio (Harvard-MGH) scanner, with further details of scan parameters detailed in both Table 1 and in prior publications (40). For Datasets 1–3, all scans were acquired using a 3T GE Discovery MR750 scanner with a 32-channel head–neck imaging coil at the Center for Cognitive and Neurobiological Imaging at Stanford. Participants were screened for MRI safety prior to any scanning procedures. In Datasets 1 and 2, MRI scans were acquired within 72 h of starting SNT (baseline) and within 72 h of completing the 5-d SNT treatment (posttreatment). In Dataset 3, each individual only received one scan (HCs). In each session, the participants underwent MRI scans consisting of structural- and rs-fMRI acquisitions. High-resolution structural images using GE’s “BRAVO” sequence (three-dimensional, T1 weighted) were acquired for the whole brain (field of view (FOV) = 256 × 256 mm matrix = 256 × 256 voxel slice thickness = 0.9 mm TR = 2,530 ms, TE = 2.98 ms, flip angle = 7°). The resting state scan consisted of a single 8-min eyes’ open resting scan. During the resting-state scans, participants were instructed to keep their eyes open and their attention focused on a central fixation point, which consisted of a black screen with a white fixation cross. The participants were also instructed to let their minds wander freely and to avoid any repetitive thoughts such as counting. Resting-state scans were collected with a 3× simultaneous multi-slice (i.e., multiband) acquisition echo planar imaging sequence: TR = 2,000ms, TE = 30ms, flip angle = 77°, slice acceleration factor = 3, FOV = 230 × 230 mm, matrix = 128 × 128 voxel, 1.8 × 1.8 mm 2 in-plane resolution, 87 contiguous axial slices parallel to the anterior commissure–posterior commissure line, yielding >1.4M voxels every 2 s. The head motion of the participants was minimized with the use of memory foam and inflatable padding. Participant alertness during the resting-state task was monitored using scanner video cameras.

Initial fMRI preprocessing followed convention as previously described using the 4dfp software suite (30, 31). Briefly, this included compensation for slice-dependent time shifts, elimination of systematic odd–even slice intensity differences due to interleaved acquisition, and rigid body correction of head movement within and across runs. Atlas transformation was achieved by composition of affine transforms connecting the fMRI volumes with the T2W and T1W structural images. Head movement correction was included with the atlas transformation in a single resampling that generated volumetric time series in 3 mm3 MNI atlas space. Additional preprocessing in preparation for latency analysis included spatial smoothing (6-mm full width at half-maximum Gaussian blur in each direction), voxel-wise removal of linear trends over each fMRI run, temporal low-pass filtering retaining frequencies below 0.1 Hz, and zero-meaning each voxel time series. Spurious variance was reduced by regression of nuisance waveforms derived from head motion correction and time series extracted from regions (of “noninterest”) in white matter and cerebrospinal fluid. Nuisance regressors included also the BOLD time series averaged over the whole brain (30, 31). Additionally, we employed frame-censoring with a frame-wise displacement (FD) threshold of 0.25 mm (63). To ensure that our results were not driven by artifactual head motion, we only included subjects with a minimum of 85% of frames remaining postcensoring. We picked the FD threshold of 0.25 mm as it allowed us to ensure no statistical difference between groups in the FD of frames retained in the analysis (dataset 1 active: 0.171 ± 0.04; dataset 1 sham: 0.177 ± 0.04; dataset 2: 0.0182 ± 0.04; dataset 3: 0.189 ± 0.04; dataset 4: 0.162 ± 0.05). We further verified that with an FD threshold of 0.25 mm, there is no statistically significant correlation between FD of retained frames and baseline MADRS scores (Pearson’s , corrected, N = 33) nor FD of retained frames and treatment response (Pearson’s , corrected, N = 23). Finally, across all subjects, we verified no statistically significant correlation between the temporal lag values seeded in the medial prefrontal cortex (mPFC) and ACC and FD values in retained frames (mPFC: Pearson’s , corrected, N = 134; ACC: Pearson’s , corrected, N = 134).

Our method for temporal delay analysis of rs-fMRI data has been extensively described in prior publications (30, 31, 33). In brief, to reduce the dimensionality of the latency analyses [number of region of interest (ROIs)], the gray matter mask was divided into 6 mm3 cubic ROIs, discounting any cubes containing fewer than 50% gray matter voxels. A TD matrix (Fig. 1) is then composed on the basis of pair-wise crosscorrelation functions between each pair of ROIs, where the peak delay is computed using parabolic interpolation. The TD matrix is then summarized through a lag projection (Fig. 1) that represents the mean delay between each cubic ROI and the rest of the brain. Seed-based latency maps are computed similarly, by computing the temporal delay between the mean time series of the seed-ROI and every other cubic ROI in the brain, to form a brain map. Matlab code implementing our temporal analyses for fMRI data is publicly available (https://github.com/ryraut/lag-code).

The neuroimaging statistical results reported in Figs. 1–5 were assessed on a cluster-wise basis using threshold–extent criteria computed by extensive permutation resampling across experimental condition (37, 64). Our approach to statistical inference, previously published (37, 65), uses the strategy, originally based on random field theory (66), of assigning significance to contiguous clusters of voxels exceeding the given threshold and extents. We compute Gaussianized t-statistic (Z-score) images, i.e., voxel-wise random effect analyses of difference contrasts, using standard formulae (67). To ensure that applied significance criteria are not affected by imprecise theoretical approximations (68), cluster-wise extent and Z-score thresholds are compiled by extensive permutation simulation of the null hypothesis (H0: no difference between groups). Thus, the frequency of clusters meeting particular threshold and the extent criteria is evaluated over 10,000 surrogate Z-score images, identical in all respects to the experimental result, except that fMRI scans are assigned to each group by random permutation. It should be noted that this approach to statistical inference is related to but not identical to controlling the family-wise error, which is formulated in terms of the largest cluster in an image (69). Here, all “significant” clusters have a controlled false positive rate of 5%. In Fig. 1, lag projection difference maps are subject to permutation, whereas in the remainder of the paper, seed-based lag difference maps are subject to permutation. In Figs. 3 and 5, P-values were computed for correlation values by transforming the Spearman’s rho value to a t-statistic with N-2 degrees of freedom, where N represents the number of subjects in the analysis. In cases where we have multiple comparisons of the same experimental question (e.g., two seeds in Figs. 2 and 3, two brain areas in Fig. 3, and two groups of control subjects in Fig. 4), we apply Bonferroni-correction to P-values.

Several different regions of interest are used to derive seed time series throughout Figs. 2–5. The left DLPFC and dorsal anterior cingulate seed maps in Figs. 2 and 3 A and B are each derived from the regions of interest identified in Fig. 1D. Note that whereas the seed-based statistical comparisons in Fig. 2 are post-hoc characterizations, the statistical comparisons in Fig. 3 A and B are conducted on an independent dataset (Dataset 2 in Table 1) from which the seeds were derived (Dataset 1 in Table 1). The left DLPFC analysis in Fig. 3C was conducted by temporal delays between the left DLPFC seed in Fig. 1D and a conjunction of the statistically significant regions shown in Figs. 2A and 3A. We chose to use a conjunction given the spatial similarity in the results, and the same mask was used in each individual subject. The ACC analysis in Fig. 3D was similarly conducted using the ACC seed in Fig. 1D and a conjunction of the statistically significant regions shown in Figs. 2B and 3B, with the same masks applied to each individual subject. The analyses in Fig. 4 were all conducted using the ACC seed in Fig. 1D. Finally, the ACC seed used in Fig. 5 also corresponds to Fig. 1D, and the delay between the ACC and its network partners computed in Fig. 5 D and E was derived using the same mask as was applied in Fig. 3D. As before, the same mask was applied to each individual subject.

Click here for additional data file.