Research ArticleNEUROSCIENCE

The physiological effects of noninvasive brain stimulation fundamentally differ across the human cortex

See allHide authors and affiliations

Science Advances  31 Jan 2020:
Vol. 6, no. 5, eaay2739
DOI: 10.1126/sciadv.aay2739

Abstract

Transcranial magnetic stimulation (TMS) is a noninvasive method to modulate brain activity and behavior in humans. Still, stimulation effects substantially vary across studies and individuals, thereby restricting the large-scale application of TMS in research or clinical settings. We revealed that low-frequency stimulation had opposite impact on the functional connectivity of sensory and cognitive brain regions. Biophysical modeling then identified a neuronal mechanism underlying these region-specific effects. Stimulation of the frontal cortex decreased local inhibition and disrupted feedforward and feedback connections. Conversely, identical stimulation increased local inhibition and enhanced forward signaling in the occipital cortex. Last, we identified functional integration as a macroscale network parameter to predict the region-specific effect of stimulation in individual subjects. In summary, we revealed how TMS modulation critically depends on the connectivity profile of target regions and propose an imaging marker to improve sensitivity of noninvasive brain stimulation for research and clinical applications.

INTRODUCTION

Transcranial magnetic stimulation (TMS) is a unique method to noninvasively modulate human brain activity. TMS allows testing causal relationships between brain activity and sensory-motor processing (1) and is capable of modulating complex cognition such as memory performance (2). Moreover, TMS has steadily evolved from a scientific tool to clinical application. Repetitive TMS (rTMS), i.e., stimulating the cortex with a sequence of electromagnetic pulses, has been used to monitor and ameliorate neurological (3, 4) and neuropsychiatric (5, 6) disorders. Only recently, the U.S. Food and Drug Administration approved rTMS as a therapeutic option for major depressive disorder (7).

Despite its undeniable impact, the replicability of rTMS effects varies substantially across individuals (8, 9). This might be due to the heterogeneous physiological effects of stimulation and a variety of available stimulation protocols. We here focus on inhibitory rTMS, i.e., the aim of decreasing neuronal excitability with either continuous 1-Hz stimulation or pulse bursts at theta (~5 Hz) frequency. The inhibitory character of both protocols has initially been established for the motor system (10) and since been generalized to inhibit other cortical regions (11, 12). Several studies, however, have reported increased brain activity with functional magnetic resonance imaging (fMRI) after 1 Hz (1315) and, more recently, after theta burst stimulation (1619). In summary, the inhibitory effect of rTMS seems not to fully generalize from motor to other functional areas of the brain. It is yet unclear why identical stimulation protocols have partly led to paradoxical effects, such as cortical excitation after inhibitory TMS. We hypothesized that the direction of modulation critically depends on the target’s connectivity profile and its underlying cellular composition.

Stimulating the cortical surface with TMS modulates a mixture of neuronal populations that use different neurotransmitters, perform different actions, and have different sensitivity to the stimulation (20). Cellular data show that rTMS modulates excitability of both γ-aminobutyric acid (GABA-) and glutamatergic neurons and thereby has different regional effects depending on the local cellular composition (21). Repeated 1-Hz stimulation particularly increased gene expression associated with synaptic plasticity (22) and GABA-producing enzymes (23), as well as GABAergic neurotransmission on the system level (24, 25). As the electromagnetic field of TMS spans several square centimeters of cortex, “identical stimulation protocols induce different early gene expression and not all brain regions respond equally to the magnetic stimulation” (26). To increase specificity and replicability of TMS, a cross-scale theory about brain stimulation is needed that takes regional heterogeneity and underlying neurophysiology into account.

To meet these demands, we combined macroscopic brain imaging with generative modeling of forward (usually excitatory) and backward (usually inhibitory) connections. fMRI identifies communication between two cortical areas via functional connectivity, a measure of temporal correlation between fMRI signals. On a more global level of interaction, graph theory methods capture a region’s functional integration into the overall brain graph (27). Others have used generative modeling to explain the endogenous brain activity underlying the fMRI signal. Particularly, spectral dynamic causal modeling (DCM) rests on a biophysically plausible model of coupled (intrinsic) neuronal fluctuations in a distributed neuronal network or graph (28).

We hypothesized that the key assumption about a frequency-dependent effect of rTMS does not generalize across the entire brain. On the basis of the heterogeneous connectivity profiles of higher cognitive and early sensory cortices (29), we systematically compared the cross-scale impact of identical 1-Hz stimulation across the human cortex. We analyzed the effect of stimulation on micro- and macroscale signaling by integrating computational modeling of cellular compartments with functional network integration on a global scale. Overall, our study revealed three major results: First, individual target identification is essential given the interindividual variability of the macroscopic brain architecture; second, identical stimulation of sensory or cognitive regions has opposite spreading effects across the cortex; third, spreading effects might be guided by a heterogeneous profile of feedback and feedforward connections of target regions.

RESULTS

Each of the 27 healthy participants underwent three counterbalanced rTMS-fMRI sessions on three different days (Fig. 1A). During each session, we identically stimulated a prefrontal (FRO), an occipital (OCC), and a temporoparietal control (CTR) region with the aim of modulating a cognitive, sensory, and functionally heterogeneous area. We measured brain activity with resting-state fMRI before (preTMS) and immediately after stimulation (postTMS). We aimed for short transition times (mean = 5.87 min, SD = 1.1 min) between the end of stimulation and postTMS and found no significant timing differences between sessions [F(2,44) = 3.24, P > 0.05, repeated-measures analysis of variance (ANOVA)]. We derived individual target regions from an online analysis of the preTMS data and loaded the coordinates into the TMS system for continuous neuronavigation during stimulation (Fig. 1B and table S1 for demographics and individual coordinates). We then applied low-frequency 1-Hz rTMS for 20 min outside the MRI scanner. Twenty-three participants (12 females; mean age = 25.74 years, SD = 3.22 years) were included in all analyses as we had to exclude two participants who did not complete all rTMS-fMRI sessions and two subjects where we could not identify target regions during the network analysis. Please find all raw imaging data and analysis scripts in the online repository of OpenNEURO (see References and Notes for the download link).

Fig. 1 Study design

(A) Each participant underwent three counterbalanced TMS-fMRI sessions on three different days. During each session, one target region (FRO, OCC, CTR; colored circles) was stimulated for 20 min with rTMS (1 Hz), and we acquired resting-state fMRI data during preTMS and postTMS. Colored overlays on the brain surface illustrate cortical template networks (31). (B) For each subject, we derived individual target spots (green spheres) within target areas from a functional network analysis of the preTMS fMRI data. (C) Brain slices with statistical parametric maps (PFWE < 0.05, corrected at cluster level, one-sample t tests) of the group average functional connectivity of each target region during preTMS calculated from the individual TMS targets. Color bars, t values.

First, we analyzed the quality of the fMRI data to allow for within- and between-subject comparisons. Per session, we identified a temporal signal-to-noise ratio [SNR(t): mean = 6.6, SD = 1.0] and framewise displacement (FD: mean = 0.13 mm, SD = 0.03 mm) in acceptable range (fig. S1) [see (30)] that did not differ between sessions [SNR(t): F(3,66) = 1.01, P > 0.05; FD: F(3,66) = 0.44, P > 0.05, repeated-measures ANOVAs]. We next validated that individually defined target areas reliably participated in the frontal (pink) and visual (violet) template networks (31). Figure 1C shows statistical parametric maps of voxels with significant functional connectivity with each of the target regions during preTMS [P < 0.05, family-wise error (FWE) corrected at cluster level, one-sample t tests]. In summary, we stimulated in each subject targets located properly within sensory (visual) and cognitive (frontal) networks, as well as a control region with heterogeneous connectivity to various networks.

Heterogeneous spreading effects for identical stimulation protocols

We next evaluated the spreading effect of rTMS for each target region (Fig. 2, fig. S2, and tables S2 and S3 for detailed results and coordinates). We found opposite effects after OCC-TMS and FRO-TMS with brain-wide increases (yellow voxels) and decreases (blue voxels) of functional connectivity, respectively. Figure 2A shows separate statistical parametric maps not only for the direct impact of stimulation, i.e., changes in functional connectivity of the stimulated target, but also for indirect effects, i.e., changes in functional connectivity of the nonstimulated target (PFWE < 0.05, corrected at cluster level, voxel-wise repeated-measures ANOVAs). We found changes neither in functional connectivity of the CTR target nor in any target region after CTR-TMS (PFWE > 0.05). In general, effects were more spatially constrained after OCC-TMS but widespread after FRO-TMS. Consistent with stimulation targets being located in the left hemisphere, changes in functional connectivity occurred more in left than right hemispheric voxels (Fig. 2A, violet bars). Stimulation effects were not spatially selective for the target region, or its associated network, but spread across several template networks (as illustrated in Fig. 2B, a summary figure of all results from Fig. 2A). Overall, identical stimulation had opposite effects on the functional connectivity of a sensory and cognitive region and spread to various functional networks.

Fig. 2 Opposite effect of TMS on brain functional connectivity

(A) Statistical parametric maps of significant changes in whole-brain functional connectivity after OCC- (top), FRO- (bottom), and CTR- (right) TMS (PFWE < 0.05, corrected at cluster level, voxel-wise repeated-measures ANOVAs). Color bar (yellow-blue), t values. OCC-TMS increased not only the functional connectivity of the OCC target (direct) but also the functional connectivity of the nonstimulated, FRO target (indirect). Conversely, FRO-TMS decreased functional connectivity to widespread cortical areas. Again, FRO-TMS had an impact not only on the functional connectivity of the stimulated, FRO target (direct) but also on the nonstimulated, OCC target (indirect). Neither of the stimulation protocols changed functional connectivity of the CTR target nor did CTR-TMS. Horizontal bar plots (violet) indicate the laterality effect of stimulation, i.e., the ratio of significant voxels in left and/or right hemispheres. (B) Summary of whole-brain functional connectivity changes as found in (A) overlaid on the cortical surface of a standard brain. Horizontal color bars indicate the ratio of significant voxels in each of the template networks as illustrated by colored outlines on the cortical surface.

Impact of stimulation on global functional integration

To investigate the stimulation effect beyond the pairwise connectivity of two regions, we next studied brain functional integration across the entire cortex. Consensus modularity analysis identified for each node the strength of local (z) and global (h) integration within a brain graph (see Fig. 3A and Methods). We consistently found three modules across all rTMS conditions, which were significantly more modular than comparable random networks on each of the three levels: individual functional connectivity matrices and individual and group coclassification matrices (P < 0.001, permutation testing; fig. S3). Figure 3B shows the group coclassification matrix of preTMS and a topological, force-directed representation of the data. The insets illustrate the location of the stimulation targets (yellow circles) and connected nodes according to coclassification. Across subjects, the FRO target had a significantly higher h than OCC target during preTMS (hFRO = 0.76 ± 013; hOCC = 0.65 ± 0.19; P = 0.01; Wilcoxon signed-rank test), which is indicative of the global integration capacity of the FRO target between the green and blue modules.

Fig. 3 Effect of stimulation on global functional integration

(A) Overview of consensus modularity analysis on the individual and group levels resulting in condition-specific parameters of local (z) and global (h) functional integration for each brain node. (B) Group average coclassification matrix (left) and its corresponding force-directed topological representation of preTMS (middle). Nodes with higher coclassification values are located closer to each other, and node size represents h. Node color indicates the modular affiliation, and abbreviations indicate assignment to template networks. Yellow circles indicate TMS targets, and insets highlight only those nodes with direct functional connectivity to the target nodes. (C and D) Group average coclassification matrices and topological representations after (C) OCC-TMS and (D) FRO-TMS with target nodes (yellow circles). Scatterplots of local (z) versus global (h) integration before (gray) and after (violet) stimulation. Note that only OCC-TMS increased global integration of the OCC target and of the entire graph. **P = 0.004, Wilcoxcon signed-rank test. Bar plot in (D) illustrates h during preTMS for all voxels that showed changes in pairwise functional connectivity. Baseline h was higher for voxels with spreading effects after FRO-TMS compared to OCC-TMS. ***P = 0.00004, Wilcoxon signed-rank test. n.s., not significant.

rTMS differentially changed the overall topology of the graph (Fig. 3, C and D). Upon visual inspection, OCC-TMS arranged the modules in a more balanced, equidistant configuration (round shape of the overall graph), indicating a relative increase of exchange between red and blue modules. FRO-TMS, however, moved green and blue modules further apart, indicating less consistent functional connectivity between the nodes. We also quantified the effect of stimulation on local and global integration values (scatter plots) between preTMS and postTMS. After OCC-TMS, global integration increased both of the OCC target (hpre = 0.65, hpost = 0.88; Fig. 3C, yellow circles) and of the entire graph (hpre = 0.80 ± 0.12, hpost = 0.83 ± 0.07, P = 0.004; Wilcoxon signed-rank test). FRO-TMS did impact on global integration neither of the FRO target nor of the entire graph (P > 0.05; Fig. 3D, yellow circles). The local integration (z) remained unaffected by any of the TMS interventions (P > 0.05). We also calculated hpre for all voxels that showed changes in pairwise functional connectivity as illustrated in Fig. 2. Across subjects, hpre was higher for all voxels that showed changes after FRO-TMS rather than for those after OCC-TMS (P = 0.00004, Wilcoxon signed-rank test; bar plot in Fig. 3D). This indicates that all regions affected by FRO-TMS were already more integrated at baseline, which might explain why stimulation broadly affected pairwise connections but was contained on the global level. In summary, stimulation of a sensory node with low global integration at baseline increased brain-wide interaction, possibly via direct connections to global integration nodes. Targeting a frontal node with initially strong global integration, however, had no impact on brain-wide integration, presumably as the stimulation effect was contained by the initially high level of dense connections.

Brain functional integration as a predictive marker for spreading effects of rTMS

We next tested the discriminatory power of global integration for the stimulation effect at the individual participant level. We calculated the difference between pre- and postTMS h for a whole-brain parcellation atlas (31) and subjected a total of N = 151 features to a random forest classifier to distinguish between OCC- and FRO-TMS sessions. Figure 4A shows the prediction accuracy for each class [OCC-TMS = 65%; FRO-TMS = 70%; 95% CI (confidence interval) = 53 to 80%], yielding an overall accuracy of 67%. Permutation testing (5000 iterations) indicated a significance level of P = 0.028 for the model. We replicated our classification result using a potentially more robust, linear classifier with less parameters. A support vector machine classifier yielded an overall accuracy of 65% (OCC-TMS = 65%; FRO-TMS = 65%; 95% CI = 51 to 78%; see fig. S4). Last, we explored the generalizability of our findings and calculated global integration parameters for the entire cortex. Figure 4B shows the regional distribution of h indicating a two-part segregation of the brain. Similar to the OCC target (violet circle), sensorimotor cortices and early integration regions have a lower h index (see also networks VIS, SM, and DAN). Frontoparietal cortices, covering higher cognitive networks such as SAL and CEN, have the highest h index, similar to our FRO target (pink circle). In summary, global integration shows a two-part distribution across the cortex, which might predict different response profiles of sensory and cognitive regions to low-frequency stimulation.

Fig. 4 Brain functional integration is a predictive marker for spreading effects of TMS

(A) Individual classification between OCC- and FRO-TMS based on h values of a whole-brain parcellation yielded an overall prediction accuracy of 67%. (Left) Confusion matrix with the prediction accuracies for every class. (Middle) receiver operating characteristic (ROC) curve and discrimination probability [area under the curve (AUC) of 0.68]. (Right) Model (green line) significantly (P = 0.027, permutation testing) deviated from the null distribution (black dashed line) after permutation testing. (B) (Left) Spatial distribution of h across the entire cortex and (right) averaged across template networks.

No effect of stimulation on the local level

Remote effects of rTMS might be simply related to local signal changes in the target region that will ultimately affect any functional connectivity measure with that region (32). We therefore analyzed three standardized imaging parameters of local fMRI signaling (Fig. 5). For each voxel in the target regions, we calculated the amplitude of low-frequency fluctuations (ALFF), the regional homogeneity (ReHo), and the SD of the signal time series. We found no significant effect of stimulation on any of the three parameters (P > 0.05, FWE corrected at the cluster level, voxel-wise repeated-measures ANOVAs). Bar plots illustrate average parameters for each region and condition. This indicates that stimulation rather affects remote functional connectivity than local signaling of a target region.

Fig. 5 No effect of TMS on local brain activity

Bar plots indicate group average (A) amplitude of ALFF, (B) ReHo of local functional connectivity, and (C) SD of the fMRI signal averaged across all voxels of each target region. PostTMS values (colored bars) after direct stimulation did not significantly differ from any other session (PFWE > 0.05 corrected at cluster level, voxel-wise ANOVA for repeated measures). Error bars represent the 95% CI of variance across subjects.

Biophysical modeling of global stimulation effects

In the final step, we tested with generative modeling how brain stimulation differentially affected feedforward and feedback connections among areas with functional connectivity changes. We used spectral DCM to characterize neuronal dynamics of local inhibitory and long-range excitatory connections in eight functional subnetworks of the template parcellation [R1 and R2, visual; R5 and R6, dorsal attention; R7 and R8, salience; R12 and R13, central executive network (31)]. Figure 6A (right) shows the group mean model of preTMS across all subjects after parametric empirical Bayes (PEB) procedures (33). We found a balanced architecture of reciprocal connections between occipital and parietal regions. Moreover, specific feedforward (green) and feedback (yellow) connections exist along an anatomical axis of sensory, parietal integration, and frontal cognitive areas. The model also estimated inhibitory self-connections (red). Figure 6B shows group differences in directional connectivity after OCC-TMS. Occipital stimulation increased local self-inhibition in the occipital cortex (R1: +0.13) and shifted the balance of directional connectivity to feedforward signaling toward parietal and frontal cortices [decreased feedback onto regions R1: −0.14 (from R13) and R6: −0.10 (from R7 and R13) and increased feedforward from R2: 0.10 (to R6)]. Conversely, FRO-TMS caused global uncoupling among all regions. Figure 6C shows decreased local self-inhibition in frontal cortex (R8: −0.10) and uncoupling of 11 of 16 directional pathways present during preTMS (dotted lines). Stimulation equally affected feedforward and feedback connections. This analysis extends the pattern of functional connectivity changes with a generative model about asymmetric cortical hierarchies in terms of feedforward and feedback connections. Stimulation of occipital regions fosters forward signaling of regions along the visual stream, while FRO-TMS leads to uncoupling of bidirectional pathways even in remote areas.

Fig. 6 Biophysical modeling of the stimulation effect on directional signaling

(A) We modeled directional signaling between regions with functional connectivity changes using spectral DCM. The model incorporates local self-inhibition (red), excitatory feedforward (green), and inhibitory feedback (yellow) signaling. Connectivity schemas illustrate directional signaling along an anterior-posterior (left-to-right) axis of the human cortex. During preTMS, the model indicated specific feedforward and feedback, as well as balanced (gray) connectivity pathways. (B) OCC-TMS strengthened feedforward signaling of occipital (R1 and R2) and parietal (R6) regions, while other connections remained in place (gray arrows). (C) FRO-TMS uncoupled 11 of 16 directional pathways (dotted lines) equally affecting feedforward and feedback connections. Violet numbers indicate changes in parameter estimates after PEB procedures with a posterior probability of >95%.

DISCUSSION

Neuromodulation with TMS has attracted much attention as it would offer an elegant way to noninvasively treat aberrant brain activity in neuropsychiatric disorders. Reports about low outcome and replicability, however, have partly compromised this method (9, 34). We here propose that a target region’s connectivity profile fundamentally affects the direction of modulation. This is a critical extension to the general assumption of cortical excitation or inhibition with specific stimulation frequencies (35) and to the variability of modulation observed during different brain states (36).

We found that identical 1-Hz stimulation had partly opposite effects on macroscopic network signaling for sensory and cognitive areas. Furthermore, generative modeling suggested that a heterogeneous cellular composition of local inhibition and remote excitation is responsible for these different response profiles. Last, we identified functional integration on the macroscale as a reliable parameter to predict a region’s response to low-frequency stimulation. Our findings provide experimental evidence for a heterogeneous signaling architecture in the human brain. Moreover, we suggest taking micro- and macroconnectivity into account when estimating the effect of regional stimulation. In summary, we provide a theoretical and practical framework to correctly target and increase sensitivity of brain stimulation in humans.

The most notable observation from our study is that 1-Hz stimulation of sensory and frontal cortices had opposite modulatory effects on that regions’ remote communication. In general, inhibitory rTMS is assumed to decrease neuronal excitability independent of the target location (10). This is unexpected as several groups have repeatedly reported increased rather than decreased activity in various cortical regions after either 1-Hz (1315) or theta burst stimulation (1619). Yet, a systematic comparison of local and global effects of stimulation to different functional areas has been missing so far. Our results suggest that the effect of modulation varies across the cortical surface and is rather determined by the extent of functional integration of a target area than by the frequency range of the stimulation protocol.

The most consistent result we identified across regions is that stimulation broadly spreads beyond the target area and associated functional networks. This remote impact is consistent with studies that used TMS to modulate resting-state functional connectivity. However, fMRI studies have usually focused on evaluating the stimulation effect within a certain functional system, such as the sensory (37, 38), motor (39), or default mode (13, 40) network. Such findings tend to convey the impression that TMS is capable of modulating a specific network in isolation. Our whole-brain approach revealed fMRI effects beyond the functional system of the target area, similar to reports about whole-brain effects using TMS–electroencephalography (EEG) (41). On the other hand, we did not find signal changes in areas below the stimulation coil with fMRI, while other imaging modalities, such as EEG and functional near-infrared spectroscopy (fNIRS), have been more sensitive to such local changes (35). Occipital stimulation spread to parietal and frontal regions along the dorsal visual stream (42), while frontal stimulation decreased functional connectivity to the salience network (covering the frontal target) and to networks across the entire cortex. This means that we identified a dichotomy of specific effects after sensory and broad effects after frontal stimulation, which reflects an established signaling hierarchy of divergence and convergence in the human cortex based on computational modeling (43) and tract-tracing studies (44). While rTMS is a promising tool to identify, and also modulate, particular functional pathways, it is important to note that spreading effects will not be confined to the network of interest alone.

We propose that functional integration on the macroscale is a suitable marker to predict specific response patterns to stimulation. Functional integration is a simple measure illustrating a region’s connectivity profile within the whole-brain graph. On the basis of the integration parameters of all cortical nodes, we successfully classified whether occipital or frontal stimulation was applied to the individual participant. Functional integration consistently distinguished among sensory and higher cognitive systems across the entire cortex. This suggests that other sensory and cognitive areas will respond similarly to stimulation as did the occipital and frontal cortex, respectively. The question then arises: How to integrate our results of pairwise functional connectivity and global integration? Occipital stimulation selectively strengthened functional connectivity with regions distributed along the dorsal visual stream. This suggests that occipital stimulation prepares the entire brain network for visual input and is in line with animal data about task-dependent activation of long-range projection neurons in sensory cortices (45). Our observation that frontal stimulation did not affect global integration is astonishing in light of the widespread decline of pairwise functional connectivity. However, such a compensatory effect has been predicted by virtual lesion modeling of brain graphs. Such studies revealed that highly integrated nodes are more resilient against massive loss of individual connections (46). In summary, rTMS effects manifest with varying, and even partly opposing, characteristics; yet, they can be consistently interpreted across different scales.

In the final approach, we created a generative model of directional signaling based on neuronal dynamics of local inhibitory and long-range excitatory connections. At baseline, spectral DCM revealed a dense network of reciprocal connections among occipital, parietal, and frontal regions in the preTMS data. These pathways are in line with major fiber tracts among occipital, parietal, and frontal lobes, such as the superior longitudinal, the occipito-frontal, and the inferior longitudinal fasciculi (47). Modeling the rTMS effects then indicated different impacts of stimulation on short-range inhibitory and long-range excitatory signaling in occipital and frontal areas. Locally, within occipital cortex, the model response to occipital stimulation was similar to a shift in the excitation/inhibition balance observed during visual processing (48). For example, Haider et al. (49) showed that local inhibition dominates excitation in amplitude and over time during awake visual processing. The model also suggested increased excitatory signaling onto lateral prefrontal cortex, both directly and via parietal cortex, similar to the signaling hierarchy along the ventral and dorsal visual stream (42). Frontal stimulation, on the other hand, led to uncoupling of various feedback and feedforward pathways reaching down to sensory areas. This is in line with macroscale data from neuroimaging about diverse connections of salience regions with both cognitive and sensory areas (50) and with tract-tracing studies about long-range feedback connections of frontal onto sensory cortices (51). Moreover, the model results converge with anatomical data about predominantly reciprocal connections among long-range pyramidal cells (52) and prefrontal cortex regions being among the top brain regions for both in- and out-degree connectivity (44). This might explain why frontal stimulation decreased widespread feedback and feedforward connections in our model. So far, frequency-dependent effects of TMS were explained by either modeling the electromagnetic effect at the site of stimulation (53) or by assuming identical modulation of distributed brain networks (54). On the path toward individualized interventions, Cocchi and Zalesky (55) stressed the need for computational models of region-specific connectivity: “Computational models to predict the local and distributed effects of TMS can help to offset the cost of empirical testing and assist with personalizing TMS protocols.” We here devise an extended, trans-scale model, anticipating both spreading effects and the direction of modulation for specific cortical targets.

Regarding the interpretation and comparability of TMS studies, it is important to remember that a variety of stimulation parameters exist. Here, we used 1-Hz stimulation at intensities defined by the motor threshold and studied the effects of stimulation during the 20 min after stimulation. Our findings are therefore specific for this combination of parameters and might not simply generalize across other inhibitory or excitatory protocols or different time windows. Still, few studies found comparable inhibitory effects between 1-Hz and theta-burst rTMS in humans (11, 12), which suggests that also theta-burst stimulation might induce heterogeneous responses across the cortex. Moreover, effects of theta burst stimulation seem to last for longer time periods. Neurotransmitter-related protein expression was elevated across several days after theta burst compared to only a few hours after 1-Hz stimulation (23). While the prolonged effects of rTMS are generally assumed to relate to synaptic plasticity, the detailed effects of stimulation on a region’s excitability via long-term depression or potentiation still remain a matter of investigation (1, 23, 26). It is, however, widely accepted that the regional level of neuronal excitability before and during stimulation has a strong impact on the size and direction of effect. Studies found that different baseline levels or priming a participant with different behavioral states strongly affects the amplitude or even direction of stimulation effects (36). As we found spatial heterogeneity to identical stimulation during an unconstrained, resting condition, this effect might also be related to different states of excitability across the cortex. It is therefore important to also compare identical stimulation protocols during different behavioral or brain states in the future.

Overall, our design allowed the systematic investigation of regional specificity of rTMS and addressed four main points of critique recently raised about the high variability of TMS results (1): We used (i) electric-field neuronavigation to identify functional target regions in each individual (53, 54) and (ii) neuroimaging to study the effects of stimulation beyond the stimulated network. (iii) We included stimulation of a control region (in contrast to sham stimulation) to test the functional specificity of our target areas. (iv) We lastly integrated macroscopic findings on the network level with a generative model to propose cellular mechanisms related to electromagnetic stimulation. Note that stimulation of a cognitive network spreads to the entire brain. The strong interconnectedness of a cognitive hub might yet compensate for the local impact and diminish the effect of stimulation. Repeated applications are therefore necessary to achieve a prolonged effect in areas with high integration capacity (7).

METHODS

Participants

Twenty-seven healthy participants (14 females; mean age = 25.56 years, SD = 3.01 years; table S1), right-handed and without any psychiatric condition, were informed of the objectives and potential risks of the study and signed a written consent inform. The study was approved by the local institutional review board and was conducted in accordance with the Declaration of Helsinki.

Brain stimulation

rTMS was delivered using an electric-field–navigated Nexstim eXimia system (version 4.3; Nexstim Plc, Helsinki, Finland) and a biphasic figure-of-eight stimulation coil. Before any rTMS session, the neuronavigation system was set up by co-registering the participants’ head to their structural MRI data [T1-weighted three-dimensional (3D)–turbo field echo (TFE) sequence], allowing to continuously track the coil position with respect to the individual target region via infrared cameras. We identified individual resting motor thresholds (rMTs) according to the maximum likelihood algorithm by mapping the cortical representation of the right abductor pollicis brevis muscle using surface muscle electrodes (Neuroline 720; Ambu, Ballerup, Denmark) and an integrated electromyography device. rTMS target regions were individually identified during a network analysis (see below) of the preTMS fMRI data and then overlaid onto structural images to guide the stimulation. Low-frequency rTMS with a frequency of 1 Hz and a stimulation intensity of 100% of the individual rMT (mean rMT = 34.4%, SD = 7.5%) was applied for 20 min to each of the target regions (1200 rTMS pulses in total) in a magnetically shielded room next to the MRI scanner. During stimulation, the coil was positioned perpendicular in relation to the skull surface, and an anterior-posterior orientation of the induced electric field was maintained using an adjustable coil holder during stimulation application.

Image acquisition

MRI data were acquired on a 3T Philips Ingenia MRI scanner using the body coil for transmission and the 32-channel head coil for signal reception (Philips Healthcare, Best, The Netherlands). We acquired multiband fMRI data during each pre- and postTMS condition (40 slices; multiband factor, MB = 2; SENSE factor, s = 2; repetition time, TR = 1250 ms; echo time, TE = 30 ms; flip angle, FA = 70°; field of view, FOV = 192 mm × 192 mm; matrix size = 64 × 64; voxel size = 3 mm × 3 mm × 3 mm). Each fMRI run lasted for 12.35 min, during which 600 functional volumes were acquired. In addition, we acquired a T1-weighted 3D-TFE during each session (170 slices; repetition time, TR = 9 ms; echo time, TE = 3.98 ms; flip angle, FA = 8°; field of view, FOV = 256 mm × 256 mm; matrix size = 256 × 256; voxel size = 1 mm × 1 mm × 1 mm). During fMRI acquisition, the scanner room was dimmed, and the participants were asked to stay awake with their eyes open. This was ensured by monitoring their eyes with an MR-compatible infrared camera (12M, MRC Systems, Heidelberg, Germany) attached to the coil.

Image data processing and analysis

We deposited the raw imaging data and analysis scripts in the online repository of OpenNEURO (https://openneuro.org/datasets/ds001927/versions/2.0.2) to allow for replication and further analyses. We performed preprocessing of the structural and functional MRI data using version 0.392 of the configurable pipeline for the analysis of connectomes (C-PAC, https://fcp-indi.github.io/). The brain representations in the figure panels were created with the help of the Nilearn library (nilearn.github.io) and the Caret software (brainvis.wustl.edu).

Preprocessing structural images

The structural images were skull-stripped using AFNI-3dSkullStrip, segmented into three tissue types using Oxford Centre for functional MRI of the brain (FMRIB) software library (FSL) for automated segmentation tool (FAST), and constrained into the individual participant tissue segmentations from standard space provided by FSL. They were then normalized to the Montreal Neurological Institute (MNI) 152 stereotactic space (2 mm isotropic) with linear and nonlinear registrations using FMRIB’s linear image registration tool (FLIRT) and FMRIB’s non-linear image registration tool (FNIRT), respectively.

Preprocessing functional images

The functional images of each run were realigned, motion-corrected to the average image using AFNI-3dvolreg, and then skull-stripped using AFNI-3dAutomask. Subsequently, the global mean intensity was normalized to 10,000, the nuisance signal was regressed, and the signal was bandpass-filtered (0.01 to 0.1 Hz). Furthermore, the preprocessed images were registered to the structural space with FSL-FLIRT using a linear transformation based on the white matter boundary information derived from the white matter tissue segmentation. The nuisance signal regression step modeled the scanner drift using quadratic and linear detrending, while the physiological noise was modeled using the five principal components with the highest variance from a decomposition of white matter and cerebrospinal fluid voxel time series (CompCor), which were derived from the prior tissue segmentations transformed from anatomical to functional space. Furthermore, head motion was modeled using the 24 regressors derived from the parameters estimated during motion realignment based on the Friston-24 Parameters, the six head motion parameters, and their 12 corresponding squared values. If not stated otherwise below, fMRI analyses were all performed in individual space. Only results were later transformed into MNI space for group statistics by applying each individual’s MNI transformation parameters of the structural image to the results maps.

Identification of individual target regions

We exported the fMRI data to an external computer, while the remaining MRI protocol of preTMS was still running. We ran an independent component analysis (ICA) of the preTMS data using FLS multivariate exploratory linear optimized decomposition into independent components (MELODIC) and decomposed the data into 17 spatial components. We then calculated the spatial cross-correlation between each individual ICA map and the following target networks from a template set of 17 networks by Yeo et al. (31): for FRO: IC_08; for OCC: IC_01; for CTR: IC_14. From each of the matching ICA maps, we then extracted the left hemispheric target areas covering the dorsoanterior prefrontal cortex (as FRO target), the occipital pole (as OCC target), and of the superior temporal gyrus (as CTR target), backprojected them into individual space, and transferred them as TMS targets onto the navigated TMS system for immediate stimulation (table S1).

Functional connectivity analysis

Functional connectivity was calculated as the pairwise Pearson’s correlation between the average time series of voxels in a seed region and the time series of all other gray matter voxels. To obtain least distorted connectivity patterns for each subject, we used individual TMS target coordinates as seed regions (5-mm-radius spheres) and calculated functional connectivity patterns in individual subject space. The preTMS functional connectivity patterns of each subject were averaged to achieve a robust pattern of baseline functional connectivity for each individual. Before we applied spatial statistics on the group level, individual functional connectivity maps were registered to the standard MNI space, z score–transformed, and spatially smoothed with a Gaussian kernel with an full width at half maximum (FWHM) of 6 mm3.

Consensus modularity analysis

Consensus modularity analysis, a graph theoretical method, identifies a unified partitioning of several graphs into nonoverlapping clusters of nodes, i.e., modules (5658). This analysis was implemented with MATLAB 2015b and the Brain Connectivity Toolbox (v2017-15-01, https://sites.google.com/site/bctnet/). We visualized the modularity decomposition using the force-directed layout representation ForceAtlas2 from the software Gephi (v0.9.2, https://gephi.org/). First, we calculated an unthresholded functional connectivity matrix on the individual level (per condition × session) using an independent parcellation atlas previously used for brain graph analysis (50). From this atlas, we selected those nodes (5-mm-radius spheres) that shared most of their voxels with both our group mask of gray matter and with any of the cognitive template networks. For each of the remaining 79 nodes, we extracted the average time series and created a cross-correlation matrix based on the pairwise Pearson’s correlation coefficient. We zeroed negative correlation values, as well as values between nodes located within a 20-mm radius (50), and applied a Fisher z-transformation. Next, we ran a consensus modularity analysis on the individual level. We iteratively (1000 times) applied the Louvain algorithm for community detection and created an individual coclassification matrix representing the frequency with which nodes were coclassified into the same module. Following recommendations in prior reports (58), we chose a threshold of τ = 0.4 for the consensus partition and iterated the process 100 times. We also present results for a range of τ values in fig. S5 and validated the consistency of our findings with more repetitions (1000 times). Last, we subjected individual coclassification matrices to a second-level consensus modularity analysis using identical parameters as above. The output of this analysis were group coclassification matrices per condition and two modularity parameters, classification consistency (z) and classification diversity (h), that illustrate a node’s local and global integration within the overall graph (56, 57).

Local signal analysis

We calculated the ALFFs, ReHo, and the SD of the signal time series for each voxel within the stimulation region. For the analysis, we used the preprocessed fMRI data with the following exceptions: The ALFF maps were calculated by computing the total power within the frequency range between 0.01 and 0.1 Hz of the nonfiltered data, whereas the ReHo maps were calculated on the basis of Kendall’s correlation between each voxel’s time series and the time series of the 27 voxels in contact with that voxel. Both measures were calculated in original space and subsequently transformed into the standard MNI space and spatially smoothed by a Gaussian kernel with an FWHM of 6 mm3.

Spectral DCM

The DCM analyses were conducted using DCM12 implemented in the SPM12 (revision 7279, www.fil.ion.ucl.ac.uk/spm/). First, we projected voxel patterns of significant functional connectivity changes back onto the template networks and identified eight subnetworks being affected by rTMS [R1 and R2, visual; R5 and R6, dorsal attention; R7 and R8, salience; R12 and R13, central executive network (31)]. Next, we extracted average fMRI signal times series for each subnetwork from individual subject space and specified a fully connected DCM model for each participant to compare all possible nested models of network interactions (33). The model was then estimated using spectral DCM, which fits the complex cross-spectral density using a power law model of endogenous neuronal fluctuations (28).

Statistical analysis

Mass-univariate voxel analysis

We performed voxel-wise group statistics applying one-sample t tests or one-way repeated-measures ANOVAs to the parameter maps of functional connectivity, ALFF, ReHo, and SD using SPM12 (www.fil.ion.ucl.ac.uk/spm/). We configured a “flexible factorial design” in SPM12 with “participants” as between-participant factor and “condition” (levels: preTMS, FRO-TMS, VIS-TMS, and CTR-TMS) as within-participant factor. Statistical testing was limited to voxels within an average gray matter mask derived from all participants. We reported significant changes with an FWE rate of P < 0.05 at the cluster level (height threshold of P = 0.001).

Modularity analysis

We validated the modular decomposition using permutation testing on three levels: (i) the individual unthresholded functional connectivity matrix and (ii) the coclassification matrix of each participant, as well as (iii) the coclassification matrix at the group level (56). On the individual level, we created random matrices matching the empirical matrices in degree, strength, and sign distribution per participant and applied the identical modularity decomposition as described above. This process was repeated 5000 times, generating a null distribution of median Q values against which we compared the magnitude of the observed sample median Q per condition (56, 59). A Wilcoxon signed-rank test (P < 0.05) was used to evaluate the effect of TMS on local and global integration parameters.

Classification

For classification, we used the random forest implementation from the scikit-learn library (v0.17.1, https://scikit-learn.org/). As features, we included the nodal difference of h between the pre- and postTMS data for all cortical nodes (N = 151; except class “undefined”) of the parcellation atlas by Power et al. (50). Critically, we used all cortical nodes of the atlas not just nodes with significant changes in our prior analysis steps, thereby avoiding any bias in feature selection (60). The calculation of h was based on the individual coclassification matrices and the power network assignments as module affiliation (50), thereby avoiding any leakage of information from the test to the training data. We evaluated the performance of the classifier using (i) a nested cross-validation (leaving out the two observations corresponding to one subject for testing) and (ii) an inner validation approach for the hyperparameter optimization of the random forest classifier (using a sequential model–based optimization implemented by the Scikit-Optimize library; skopt, https://github.com/scikit-optimize/scikit-optimize), iteratively tuning the following parameters following the recommendations by Probst et al. (61): maximum depth of the tree, number of features, minimum number of samples, and minimum number of samples required to be at a leaf node. We statistically validated the observed accuracy using permutation testing (P < 0.05, 5000 iterations) randomizing the class labels.

PEB framework for DCM

The subject-specific DCMs were taken to the second level where we used PEB routines for group level inference (33); these routines assess how individual (within-subject) connections relate to group means, taking into account both the expected strength of each connection and the associated uncertainty. Specifically, we created three separate second-level PEB models to examine directional connectivity at baseline (preTMS) and changes after OCC-TMS and FRO-TMS within eight functional subnetworks. Next, we used Bayesian model reduction to test all nested models within each full PEB model [assuming that a different combination of connections could exist for each participant (33)] and to “prune” connection parameters. The parameters of the best 256 pruned models were averaged and weighted by their evidence (Bayesian model averaging) to generate group estimates of connection parameters. Last, we compared models using free energy and calculated the posterior probability for each model as a softmax function of the log Bayes factor. We report effects (connection strengths) as significant with a posterior probability of >0.95.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/5/eaay2739/DC1

Table S1. Participants demographic information and individual TMS target coordinates in MNI space.

Table S2. Statistical significant changes in whole-brain functional connectivity after OCC-TMS as presented in Fig. 2 and fig. S2A.

Table S3. Statistical significant changes in whole-brain functional connectivity after FRO-TMS as presented in Fig. 2 and fig. S2B.

Fig. S1. Image quality assessment.

Fig. S2. Cross-sectional representation of Fig. 2A.

Fig. S3. Statistical significance of the modularity results.

Fig. S4. Classification results using linear support vector machine (SVM).

Fig. S5. Effect of the parameter τ on the global functional integration results.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: Funding: This study was conducted with financial support to V.R. from the Deutsche Forschungsgemeinschaft (DFG, grant agreement no.: 273427765) and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no.: 759659). G.C. received a scholarship from the Colombian Government (COLCIENCIAS scholarship, grant agreement no.: 646, 2014). A.R. is funded by the Australian Research Council Discovery Early Career Research Award Fellowship (DE170100128) and the Wellcome Trust. Author contributions: Conceptualization: V.R. Methodology: G.C. and A.R. Formal analysis: G.C. Investigation: G.C., N.S., and K.K. Resources: S.M.K. Data curation: G.C. and V.R. Writing (original draft): G.C. and V.R. Writing (review and editing): N.S., A.R., and S.M.K. Visualization: G.C. and V.R. Supervision: V.R. Competing interests: N.S. received honoraria from Nexstim Plc (Helsinki, Finland). S.M.K. is a consultant for Brainlab AG (Munich, Germany), Spineart Deutschland GmbH (Frankfurt, Germany), and Nexstim Plc (Helsinki, Finland) and received honoraria from Medtronic (Meerbusch, Germany) and Carl Zeiss Meditec (Oberkochen, Germany). The authors declare no other competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Imaging data and analysis scripts can be downloaded from the online repository OpenNEURO (dataset DOI: 10.18112/openneuro.ds001927.v2.0.2, https://openneuro.org/datasets/ds001927/versions/2.0.2) to allow for replication and further analyses.
View Abstract

Stay Connected to Science Advances

Navigate This Article