Modulation of Long-Range Connectivity Patterns via Frequency-Specific Stimulation of Human Cortex.

There is increasing interest in how the phase of local oscillatory activity within a brain area determines the long-range functional connectivity of that area. For example, increasing convergent evidence from a range of methodologies suggests that beta (20 Hz) oscillations may play a vital role in the function of the motor system [1-5]. The "communication through coherence" hypothesis posits that the precise phase of coherent oscillations in network nodes is a determinant of successful communication between them [6, 7]. Here we set out to determine whether oscillatory activity in the beta band serves to support this theory within the cortical motor network in vivo. We combined non-invasive transcranial alternating-current stimulation (tACS) [8-12] with resting-state functional MRI (fMRI) [13] to follow both changes in local activity and long-range connectivity, determined by inter-areal blood-oxygen-level-dependent (BOLD) signal correlation, as a proxy for communication in the human cortex. Twelve healthy subjects participated in three fMRI scans with 20 Hz, 5 Hz, or sham tACS applied separately on each scan. Transcranial magnetic stimulation (TMS) at beta frequency has previously been shown to increase local activity in the beta band [14] and to modulate long-range connectivity within the default mode network [15]. We demonstrated that beta-frequency tACS significantly changed the connectivity pattern of the stimulated primary motor cortex (M1), without changing overall local activity or network connectivity. This finding is supported by a simple phase-precession model, which demonstrates the plausibility of the results and provides emergent predictions that are consistent with our empirical findings. These findings therefore inform our understanding of how local oscillatory activity may underpin network connectivity.

Results

20 Hz tACS Does Not Alter Resting Activity in the Stimulated M1

We first wished to investigate whether 20 Hz transcranial alternating-current stimulation (tACS) would induce changes in cortical activity (Figure 1). A whole-brain voxel-wise analysis indicated that, as predicted, there was no net change in activity during 20 Hz stimulation compared with either 5 Hz tACS or sham tACS (thresholded at Z > 2.3, p = 0.05 (corrected)). Additionally, there were no activity changes during 20 Hz tACS within our pre-specified regions of interest (ROIs; left primary motor cortex [M1], right M1, left premotor cortex [PMC], and right PMC; see STAR Methods for details; repeated-measures ANOVA [RM-ANOVA]: no main effect of ROI [F(2,22) = 3.093, p = 0.071], no main effect of tACS condition [F(2,22) = 1.920, p = 0.17], and no interaction between ROI and stimulus condition [F(6, 66) = 0.907 p = 0.495]).

Figure 1

Experimental Outline

Each participant had three scans, acquired on the same day, and during which 20 Hz, 5 Hz, or sham tACS was applied with a sinusoidal waveform and no current offset to the left M1, with the order counterbalanced across the group. Real tACS was performed for 60 s, with 10 s ramp-up and ramp-down periods on either side, repeated four times in each run. Sham stimulation consisted of 10 s ramp up and ramp down only. There were 110 s between each stimulation period. Subjects were advised to keep their eyes open and to look at a cross displayed centrally throughout the scans. Insets show details of the current amplitude for each stimulation period. Participants were asked to rate levels of paraesthesia, pain, and phosphenes on visual analog scales between each scan. See also Figure S1.

20 Hz tACS Does Not Modulate Overall Motor Network Connectivity

We primarily wished to investigate the effects of M1 tACS on functional connectivity across the motor network. To assess the strength of connections between key network nodes, we used both ROI seed-based analyses and an independent-component analysis (ICA) approach, which identifies spatially distinct networks of regions with correlated resting blood-oxygen-level-dependent (BOLD) activity, the so-called resting-state networks (RSNs) [16, 17, 18]. For the ICA analyses, we were primarily interested in the motor network and additionally identified the default mode network (DMN), a well characterized network that does not include M1, as an anatomical control. Network strength gives a measure of functional connectivity within that network and has been shown to be a sensitive metric for better communication between the major network nodes [19, 20].

We first performed an ROI analysis to investigate changes in functional connectivity during tACS between the stimulated (left) M1 and other major network nodes (right M1; left PMC). There were no significant changes in M1 connectivity with 20 Hz tACS compared with the other frequencies (RM-ANOVA, one factor of frequency [20 Hz, 5 Hz, sham] and one factor of ROI [right M1, left PMC]: main effect of frequency, F(2,22) = 0.87, p = 0.431; frequency × ROI interaction, F(2,22) = 8.13, p = 0.002; M1-M1 connectivity shown in Figure 2A). We went on to explore the frequency by ROI interaction and demonstrated no frequency-specific changes in M1-M1 connectivity (RM-ANOVA, one factor of frequency [20 Hz, 5 Hz, Sham], F(2,22) = 2.30, p = 0.13) but a significant decrease in M1-PMC connectivity during 20 Hz tACS (RM-ANOVA, F(2,22) = 9.51, p = 0.001; post hoc t test [20 Hz v Sham], t(11) = 4.1, p = 0.002; Figure S2A).

Figure 2

Modulation of the Relationship between M1-M1 Connectivity and Motor Network Connectivity

(A) ROI analysis showed 20 Hz stimulation did not change M1-M1 connectivity compared to either 5 Hz or sham stimulation. Bars indicate mean ± SEM.

(B and C) Similarly, there was no change in overall network strength in either the motor network (B) or default mode network (DMN) (C) as a result of stimulation. Bars indicate mean ± SEM.

(D–F) Relationship between motor network strength and M1-M1 connectivity. The expected close relationship between M1-M1 connectivity and motor network strength, seen with both sham (F) and 5 Hz (E) stimulation was lost with 20 Hz stimulation (D), suggesting that the pattern of connectivity within the motor network was significantly changed by local stimulation at the beta frequency. Results of linear regression and 95% confidence limits shown in (E) and (F).

(G–I) As expected, the relationship between M1-M1 connectivity and RSN strength was anatomically specific, with no correlation between M1-M1 connectivity and DMN strength. As Pearson’s correlation coefficient is not normally distributed, an r-to-Z transformation was performed for all measures of M1-M1 connectivity. Network strength is calculated as the mean parameter estimate across the network and is given in arbitrary units.

The asterisk indicates significant difference (p < 0.05) between sham stimulation and 20 Hz stimulation in the relationship between M1-M1 connectivity and motor network strength. See also Figure S2.

We next investigated tACS-related changes in connectivity within subject-specific motor network and DMN maps (Figure 3; see STAR Methods for details of subject-specific networks derivation [21]). There were no frequency-specific changes in network strength in the motor network or the DMN (RM-ANOVA, one factor of frequency (20 Hz, 5 Hz, sham), one factor of network (motor, DMN); main effect of frequency, F(2,22) = 2.902, p = 0.76; frequency × network interaction, F(2,22) = 2.083, p = 0.148; Figures 2B and 2C).

Figure 3

Resting State Networks

ICA-derived group mean (A) motor resting state network and (B) default mode network.

20 Hz tACS Uncouples Components of the Motor Network

It has previously been demonstrated that M1-M1 connectivity and motor RSN strength are ordinarily closely related [21]. Our primary hypothesis was that stimulation of the left M1 at its resting resonant frequency would entrain left M1, together with its strongest reciprocal connections, to the extrinsic rhythm. The inter-areal connectivity, as judged by functional MRI (fMRI), would not necessarily be expected to change, but the phase of activity would be drawn away from the steady-state phase relationships established in the remainder of the motor RSN. Therefore, the strong inter-areal connectivity between the M1s would be uncoupled from the remainder of the motor RSN.

As expected, there was a strong correlation between M1-M1 connectivity and motor network strength in the sham condition (R2 = 0.74, p < 0.001; Figure 2F). However, this relationship was lost during 20 Hz stimulation (R2 = 0.05, p = 0.47; 20 Hz versus sham, Fisher’s r-to-Z: Z = 2.25, p = 0.024; Figure 2D). Critically, this uncoupling was frequency selective and not observed with 5 Hz tACS (5 Hz: R2 = 0.74, p < 0.001; 20 Hz versus 5 Hz, Fisher’s r-to-Z: Z = 2.23, p = 0.026; Figure 2E). No significant correlation was observed between M1-M1 connectivity and DMN strength at any tACS frequency (Figures 2G–2I).

To confirm an uncoupling between left M1 and the rest of the motor network, we additionally performed an ROI analysis to investigate connectivity between the left M1 and the motor network excluding both M1s (rest of network). This confirmed a frequency-specific decrease in connectivity between left M1 and the rest of the network during 20 Hz stimulation compared with sham (Figure S2B, RM-ANOVA, main effect of stimulation: F(2,22) = 4.243, p = 0.02; follow-up t test between 20 Hz and sham: t(11) = 3.168, p = 0.009).

However, entraining M1 oscillations with 20 Hz tACS might be expected to act as noise in circuits that are not strongly connected to M1 and may not mirror its resonance profile. We therefore calculated the strength of the connectivity in areas within the motor network but outside the M1s. There was a significant decrease in connectivity in these areas with 20 Hz compared to sham (t(11) = 2.31, p = 0.04), supporting the hypothesis that the tACS-entrained oscillations in M1 were insufficient to fully capture phase relationships within secondary motor areas but still managed to weaken existing phase alignments favoring communication.

Finally, we wished to investigate whether there were any systematic differences in signal to noise ratio (SNR) across our three scan sessions that might explain our results. We therefore calculated the SNR in each scan session for each subject separately and demonstrated no significant differences in SNR between any of the three stimulation conditions (thresholded at Z > 2.3, p = 0.05 [corrected]).

A Simple Phase-Precession Model Explains Key Features of Our Data

Taken together, our findings point toward a system of coupled brain regions that rely on temporal channels of synchronous bursting for communication, and that moderately weak stimulation (tACS) close to the resonant frequency of this system acts to actively decouple the constituent elements of the network. We wanted to demonstrate the feasibility of such a scenario by constructing a mathematical model that relies only on the principles of phase precession (Figure 4A).

Figure 4

Phase-Precession Model

(A) Model topology.

(B–F) The model was able to faithfully reproduce key features of the observed experimental data, including (B) the (frequency-specific) loss of association between M1-M1 and overall motor network strength and (C and D) the reduction in net coupling of the secondary motor system, despite relatively preserved M1-M1 coupling. Our model additionally predicted that a M1-M1 loop with intrinsically weak connectivity is more readily entrained by 20 Hz (but not 5 Hz) tACS, whereas a more synchronized loop is less affected and may even begin to desynchronize (E). Baseline M1-M1 connectivity is plotted on the x axis, change in M1-M1 connectivity due to 20 Hz tACS on the y axis. This was confirmed by a secondary analysis of the fMRI data (F). Δ, change in PSI due to stimulation.

This model was able to faithfully reproduce key features of our experimental data. These included: (1) the frequency-specific loss of association between M1-M1 and overall motor network strength (R2 = 0.913 without stimulation versus R2 = 0.130 during 20 Hz stimulation; difference of correlation test, z = −10.5, p < 0.001; 5 Hz stimulation had no discernible effect, R2 = 0.902, z = −0.430, p = 0.667); (2) a strong reduction in net coupling of the secondary motor system (Cohen’s d = −5.183), absent with 5 Hz stimulation; and (3) relatively preserved M1-M1 coupling (Cohen’s d = 1.391; Figures 4B–4D).

Our model additionally predicted that a M1-M1 loop with intrinsically weak connectivity is more readily entrained by 20 Hz (but not 5 Hz) tACS than a more synchronized loop, which may even begin to desynchronize (Figure 4E). This was found to be the case in a secondary analysis of the fMRI data, where we related M1-M1 connectivity during sham stimulation, as a baseline metric, to the change in M1-M1 connectivity between 20 Hz tACS and sham stimulation (R2 = 0.52, p = 0.008; Figure 4F). Note that while this model captures (even predicts) many of the characteristics observed in our experimental data, it represents a single topology tuned independently of the data itself and employs phase synchronization as a proxy for correlations in BOLD activity. As such, the dynamic range of the model will be far narrower than that expected from the experimental data, which is considered across subjects.

Discussion

This study was performed to investigate the role of M1 beta oscillations in determining the network-level functional connectivity of the region. tACS is believed to superimpose alternating subthreshold windows of excitation and inhibition on underlying neurons, with this effect being more pronounced through local resonance when stimulation is delivered at the relevant frequencies [22]. Accordingly, we delivered tACS at 20 Hz over the M1 at rest with the intention of entraining local beta activity, thereby shifting the phase of beta oscillations away from that dictated by the steady-state dynamics of the wider motor network. Although weak currents would not be predicted to change the overall level of local activity, they might be proposed to draw M1 and its strongest reciprocal connections away from those phase relationships with the remainder of the motor network that would otherwise be optimal for communication [23]. This effect might by induced by synchronization of the neural firing to the extra-cranially applied sinusoidal current as has been shown in a rat model [24], and increases in beta activity have been demonstrated after 20 Hz rTMS [14].

We tested 20 Hz, 5 Hz, and sham tACS applied to the left M1 during resting state fMRI and found that while beta-frequency tACS led to no significant overall change in M1-M1 connectivity or in motor network strength, it changed the relationship between these two metrics. Previous work has established a significant relationship between connectivity between the left and right M1 and connectivity within the motor network as a whole [21, 25], reflecting the role of the M1s as the major nodes within this network. Although the study was limited to 12 participants and therefore might not have detected subtle changes in connectivity, the finding that beta-frequency tACS significantly modulates this relationship is important, as it suggests that changing the phase of local activity at the resonant frequency of M1 can modulate the connectivity pattern of that area. It should be noted that a direct effect of the reference electrode on the supra-orbital ridge cannot be ruled out. Similarly, it is possible that some of the current from the electrode positioned over M1 spread to surrounding cortical regions, such as the premotor cortices. However, the M1-specific findings demonstrated here are consistent with the major effects of stimulation being maximal over this site, perhaps because the 20 Hz stimulation frequency was deliberately chosen to be close to the local resonance frequency.

A parsimonious explanation for the results seen here might be that the substantial structural connections between the two M1s and their similar resonance characteristics means that the phase of both left and right M1s becomes relatively entrained by 20 Hz tACS, leading to greater maintenance of phase coupling and therefore correlated fluctuations in BOLD signal between these primary motor areas. However, the capture of phase by extrinsic tACS uncouples M1-M1 connectivity from the remainder of the motor network so that the fluctuations in BOLD signal between the two subcomponents are no longer correlated. In practice, though, it is likely that MI activity also influences phase relationships, and hence communication, between non-primary motor regions. Hence, on closer inspection, when connectivity within secondary motor areas was considered alone, uncoupling of M1 from the remainder of the motor network did decrease network connectivity in other motor regions. Presumably, in this case, connections with M1 are weaker and resonance characteristics less well matched, so that M1 entrainment acts as a source of noise, partially disrupting those optimal phase relationships for communication that may spontaneously emerge within secondary motor areas [26]. The finding that 20 Hz tACS decreases M1-PMd connectivity in a manner that significantly co-varied with the effects of 20 Hz tACS on M1-M1 connectivity further supports the conclusion that tACS at 20 Hz entrained the phase of underlying M1 oscillations, which in turn affects the phase relationships between connected regions. In order to explore the plausibility of this hypothesis, we developed a simple-phase-precession model, which was able to explain key aspects of our data, as well as provide emergent properties that we were able to support with our empirical data.

It is important to note that the functional connectivity modulations were quick to emerge and observed during relatively short periods of tACS. We specifically designed this experiment to avoid tACS aftereffects as we wished to focus on entrainment effects during stimulation, and we use repeated stimulation periods within the same scanning period to reduce across-session variation in BOLD fluctuations. Convincing tACS-related behavioral effects have previously been observed during short-periods of tACS [9, 11], leading to the speculation that these effects may be underpinned, at least in part, by the functional connectivity changes observed here.

Communication through Coherence Hypothesis

We directly tested a commonly held theory that proposes that the efficacy of information flow between different brain areas depends on the phase alignment between their intrinsic activities, the so-called “communication through coherence” hypothesis. According to this hypothesis, optimal phase alignment between functionally connected cortical areas predicts an increase in correlation between the two areas, where correlation is a proxy for information exchange [6, 7, 27]. In practice, such optimal phase alignment may be the natural result of the steady-state dynamics of interconnected neural circuits [26]. Our findings, while necessarily indirect to avoid confounds that can occur with electroencephalographic (EEG) and magnetoencephalographic (MEG) signals due to electrical stimulation [12, 25], nevertheless provide empirical support for the “communication through coherence” hypothesis. Further, our data suggest that the impact of phase capture of one area upon phase relationships in interconnected regions may be dependent on the strength of recurrent connectivity and resonance characteristics. It will be important for future studies to confirm these findings with EEG or MEG approaches and to determine the extent to which task-dependent changes in recurrent connectivity and local resonances modify entrainment.

Other Applicable Hypotheses

The data presented here are also open to interpretation in the context of other theories of inter-areal communication. For example, it has been proposed that power-power coupling between areas is crucial for communication [26]. In that MEG study power envelope correlations occurred on the timescale of several seconds. These therefore overlap with the frequencies of BOLD covariation detected in the present study, raising the intriguing possibility that the previously reported power-power coupling and the present findings may be related. Indeed, if one assumes that phase synchronization across neural populations leads to increased post-synaptic efficacy through temporal integration and thereby improved communication, then low-frequency dynamics in phase synchronization should be translated in to low-frequency amplitude covariations. Whether phase synchronization (coherence) or power coupling are the mechanistically important aspect supporting neural communication then becomes an issue for empirical testing [27].

The Role of Beta Oscillations

The beta band frequency forms the dominant oscillatory activity in the primary motor cortex at rest [3, 4]. In general, inter-areal connectivity, quantified here via resting-state fMRI, is thought to be primarily driven by alpha and beta oscillatory activity within key network nodes [28, 29, 30, 31, 32], evidence supported by a recent study demonstrating that alpha-frequency tACS led to significant modulation of functional connectivity [33], though necessarily this evidence is indirect.

However, this does not exclude a role for activity at other frequencies, particularly the gamma band, in network-level connectivity [28, 34]. Indeed, it might be that stimulating at the gamma frequency would lead to a change in the overall connectivity within the network, especially during movement, and this should be investigated directly in future studies.

This study is, to our knowledge, the first to explore the role of modulating M1 beta phase in network-level connectivity and lends substantial support to the hypothesis that phase coupling between anatomically distant regions underpins functional connectivity and therefore communication, between them.

STAR★Methods

Key Resources Table

Contact for Reagent and Resource Sharing

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Dr Charlotte Stagg (charlotte.stagg@ndcn.ox.ac.uk).

Experimental Model and Subject Details

12 healthy volunteers (7 male; mean age: 26 (range 21-42)) gave their informed consent to participated in the study in accordance with local ethics committee approval (NRES Committee South Central - Berkshire B; 13/SC/0413). Subjects were right-handed [35], had no previous neurological or psychiatric history and no contraindications to transcranial stimulation. Volunteers participated in two MRI sessions. The first session comprised three resting state scans, each lasting approximately 14 min. During each scan, four tACS periods of the same condition, each lasting 80 s (including 10 s ramp-up and ramp-down) were applied, with 110 s tACS-free periods between them (Figure 1). Current was switched off between the ramps for sham stimulation. The order of scans was counterbalanced across subjects. Subjects were asked to keep their eyes open and fixate on a cross presented in the center of the screen. Alertness, and possible tACS side-effects were assessed throughout the session (see Figure S1).

In order to identify functional ROIs, a second MRI session was performed.

Method Details

MR acquisition

Whole-brain functional resting-state MRI (rs-fMRI) was performed using a multi-band echo planar imaging (MB-EPI) sequence on a 3T MRI system (Magnetom Verio 3T, Siemens) using a 32-channel head coil. For registration purposes a high-resolution T1-weighted 3D structural image was also acquired.

Resting state imaging

Whole-brain functional resting-state MRI (rs-fMRI) was performed using a multi-band echo planar imaging (MB-EPI) sequence (TR: 1300 ms, TE: 40 ms, 72 × 2 mm axial slices providing whole-brain coverage, FOV: 212 mm x 212 mm, voxel size: 2x2x2 mm) [36, 37]. During acquisition, subjects were instructed to lie awake and still in the scanner while fixating on a central fixation cross. For registration purposes a high-resolution T1-weighted 3D structural image (TR: 2040 ms, TE: 4.7 ms, flip angle: 8°, field of view: 192 mm x 192 mm, voxel size: 1x1x1 mm) was also acquired for every participant.

Functional Localizer scan

This session consisted of a sequential externally guided unimanual finger tapping task. Subjects were requested to tap either with their right or left index finger on a button box at 1 Hz. The speed of the tapping was indicated via a green blinking dot. The session consisted of 6 tapping blocks interspersed with 6 rest blocks, each lasting 30 s (total scanning time approximately 6 min). Before each tapping block subjects were instructed to start tapping with either the right or left hand via a short instruction displayed on the screen.

Transcranial alternating current stimulation (tACS)

An MRI-compatible DC-Stimulator-plus (Neuroconn GmbH, Germany) was used to deliver a 1 mA current via two 5x7 cm conductive rubber electrodes. One electrode was placed over the left primary motor cortex [M1; C3 position according to the International 10-20 EEG system]; the other electrode was fixed over the right supraorbital ridge. We used a sinusoidal waveform without DC offset, resulting in a mean current density (peak-to-peak) of 0.029 mA/cm2 under the electrodes. While we did not continuously record impedance throughout the study, it was kept between 10 and 20 kΩ at all times, in line with local operating practice.

In order to assess the subjects’ alertness during the scans, the question “Are you tired”? was displayed on the screen for 5 s, 10 s before each stimulation block. Subjects were asked to respond via a hand-held button box using the index finger (“yes”) or middle finger (“no”) of their right hand in order to give a response. After each scan participants were asked to rate their experience of paraesthesia, discomfort, and the strength of phosphenes from 0 – 10 on visual analog scales with anchors (0: “no paraesthesia,” “no pain,” “no phosphenes”; 10: “high paraesthesia,” “high pain,” “strong phosphenes”; see Figure S1).

Image analysis

All data were analyzed using FSL tools (http://www.fmrib.ox.ac.uk/fsl) [38, 39]. Individual preprocessing steps included motion correction using MCFLIRT [40], non-brain tissue removal using BET [41], non-linear high-pass temporal filtering equivalent to 150 s (0.007 Hz) and spatial smoothing using a Gaussian kernel of FWHM 5.0mm. FMRI volumes were registered to T1-weighted high-resolution individual subject anatomical images using linear registration, and then to standard MNI152 template image using FMRIB’s Nonlinear Image Registration Tool (FNIRT) [39].

Next, denoising was performed using Independent Component Analysis (ICA), as implemented in FSL. Independent Component Analysis (ICA), as implemented in FSL was run for each subject and each scan separately. Each component from the ICA was classified as either signal or noise in two ways to ensure reproducibility. First, identification of artifactual components was performed manually by a blinded researcher. However, the decision of whether a component is classified as signal or as noise requires a high level of expertise of signal and noise fluctuations’ spatiotemporal characteristics, and can therefore be confounded by subjective judgment, several “full-automatic” approaches to ICA classification have been developed [42]. To ensure the reliability of our ICA classification we applied FIX (FMRIB’s ICA-based Xnoiseifier) [43]. A comparison between the automatic and manual ICA classification revealed an inter-rater co-efficient of 0.95, demonstrating the accuracy of our manual ICA classification. The artifactual components identified from our manual approach were then removed. As the model-free ICA approach does not allow the inclusion of regressors of no interest for the visual display and the subsequent button press response we excluded the volumes where these events occurred [44].

Analysis of brain activity changes

Voxel-wise analyses. Region of Interest (ROI) masks were defined as follows:

Primary motor cortex (M1): We determined the group mean activity within the left M1 during tapping of the right finger, which was then masked by a generous anatomical mask of the left M1 [45]. We then extracted the peak voxel and created a 10x10x10 mm ROI centered on these co-ordinates. The mask for the right M1 was defined similarly from the left tapping data.

Premotor cortex (PMC): We derived PMC masks from a connectivity-based parcellation of the premotor cortex [46], mapped onto the standard space MNI template.

To investigate changes in brain activity induced by tACS we used a General Linear Model (as implemented via FEAT in FSL [38]). The first level model consisted of one boxcar regressor modeling the tACS stimulation.

A second level mixed-effects analysis was used to calculate group mean maps of areas changed activity during stimulation versus baseline across all subjects for each stimulation condition separately. In addition, a second level mixed-effects analysis was performed to investigate differences in activation between stimulation conditions.

Region of interest (ROI) analyses. Given our strong a priori hypothesis as to the anatomical distribution of tACS-induced brain activity changes based on the known connectivity of M1, we also performed a ROI analysis to investigate potential effects of tACS on brain activity in the motor network.

The mean percentage BOLD signal change across all activated voxels within each ROI was calculated for each subject. BOLD changes were then compared between conditions using a repeated-measures ANOVA with one factor of tACS frequency (20 Hz, 5 Hz or sham) and one factor of ROI (left M1, right M1, left PMd and right PMd). When sphericity was violated appropriate corrections were performed.

Functional Connectivity analysis

Group ICA. In order to investigate whether tACS induced connectivity changes within brain networks the concatenated fMRI dataset was analyzed using ICA as implemented in MELODIC [17]. Data were decomposed into 20 components and RSNs of interest were identified using spatial correlations against previously defined maps [17]. A dual-regression approach was used to identify subject-specific RSN maps for each tACS condition [47, 48].

The subject-specific motor network and DMN RSN was then masked by the corresponding group mean RSN and the mean value within this region extracted for each subject, giving a measure of the strength of functional connectivity within each of the selected RSNs [21].

To assess network strength of the non-primary motor regions, the subject-specific motor network was masked first by both the group mean RSN as an inclusion mask, and then by an anatomical M1 mask as an exclusion mask [45]. The mean strength within the remaining non-primary motor regions was then calculated as above.

Seed-based correlations. Because we were primarily interested in connectivity changes between the stimulated M1 and other key nodes within the motor network we investigated connectivity changes between the left M1 and right M1 and between the left M1 and left PMd. We extracted the individual time course of the fMRI signal from the ROIs as described above. A Pearson’s correlation analysis was performed on the time-courses from the left and right M1; and the left M1 and the left PMd. As this is non-normally distributed we performed a Fisher’s r-to-Z transformation prior to running statistical analyses.

Phase-precession model

There is a growing body of evidence that implicates phase-precession as a physiologically plausible mechanism for the emergence of slow-frequency fluctuations in MEG and BOLD fMRI [33, 34]. These dynamics have been modeled using Kuramoto-style oscillators [49] that have faithfully reproduced amplitude fluctuations as emerging from metastable oscillatory dynamics. In order to demonstrate that our hypothesis is compatible with such a phase-precession framework, we employed a Kuramoto-style model with dynamics given by,where represents the phase of the ith oscillator, is the natural frequency of that oscillator, is the strength of influence of node j on node i, and represents a noise source. Our network incorporated five oscillator nodes, whose topography is shown in Figure 4A. All nodes received independent noise of standard deviation , excepting the tACS node where .

Model parameters were selected to simultaneously 1) maximize the correlation between M1-M1 and the broader motor network in the absence of stimulation, 2) minimize the correlation between M1-M1 and the broader motor network during 20 Hz stimulation, 3) reduce SM↔SM coupling during 20 Hz stimulation and 4) minimize the change in M1-M1 coupling due to 20 Hz stimulation. We are interested in the coupling between nodes reflected in their relative phase-alignments, which will be adopted as a reasonable surrogate of the low-frequency covariation in BOLD within this simulation. This is based on the premise that phase synchronization across neural populations (or coupled oscillators as here) leads to increased post-synaptic efficacy through temporal integration and thereby improved communication, so that low frequency dynamics in phase synchronization should be translated to low frequency amplitude covariations detectable in the BOLD signal [49]. To be explicit, BOLD was not used here to track beta oscillations per se. Rather we used it to follow emergent dynamics in neural activity at low frequencies arising through intrinsic modulations in the phase synchronization at higher frequencies, or due to small frequency mismatches between intrinsic beta activities and the 20 Hz drive provided by tACS. Both lead to low frequency modulations in the envelope of phase synchronization. In the case of tACS this modulation will have a frequency that is related to the mismatch between the frequency of intrinsic beta activities and the beta tACS. The stimulation frequency at 20 Hz was carefully chosen so that it was very close to the peak frequency of beta (mean 20.1 Hz, SD = 2.07, n = 20) recorded over M1 in healthy control subjects of similar age in our lab [50], thereby constraining phase drifts to the frequency band that can be detected in the BOLD signal.

The objective of this simulation was to demonstrate, in line with our hypothesis, that the phase-precession model is capable of reproducing the behavioral phenomena observed in this study. Thus, although the model was optimized as a deterministic system (using a fixed random seed), the model nevertheless illustrates that topologies exist in a phase-precession framework that are compatible with our experimental observations and hypotheses. In addition, different random seeds were selected to generate the realizations depicted in Figure 4, thus the model appears broadly generalizable. Parameters include the natural frequencies (allowable range 18.0-22.0 Hz; except tACS which was fixed at 20.0 Hz), connection weights M1-M1 (symmetric), M1→SM, SM→M1, tACS→M1, and the standard deviation of noise . Connection weights were constrained between [-1, 0] which facilitated attractor dynamics. Optimization was conducted in two passes first through a custom written swarm optimization routine, then an interior-point algorithm. Once optimized, 100 sample simulations of 10 s were generated under fixed initial conditions but with 10 s burn-in. We quantified the net change in phase synchronization index (PSI) due to stimulation versus no stimulation using Cohen’s d, correlations by Pearson’s Product-Moment method, and differences in correlations due to stimulation by normalizing Pearson’s correlation coefficients by Fisher’s z-transform and computing the z-score of differences with associated two-tailed statistics. The model was constructed, simulated and analyzed using MATLAB (R2016a, Mathworks, USA).

The derived parameters of the model were – Frequencies: 19.7 Hz (M1(L)), 18.9 Hz (M1(R)), 20.5 Hz (SM), 21.3 Hz (SM); Connection strengths: −0.289 (M1(L) ↔ M1(R)), −0.473 (M1 → SM), −0.604 (SM → M1) and −0.500 (tACS → M1(L)); Noise contribution, .

Quantification and Statistical Analysis

Whole Brain Analyses

For all group-level analyses, Z-statistic images were thresholded using clusters determined by Z > 2.3 and a (corrected) cluster significance threshold of p = 0.05 within FSL.

Region of Interest Analyses

Statistical analyses were performed using SPSS (version 20.0) and the statistics toolbox of MATLAB (version 7.11.0.584). Connectivity changes were compared between conditions using repeated-measures ANOVA with one factor of tACS frequency (20 Hz, 5 Hz or sham) and one factor of ROI (left M1, right M1, left PMd and right PMd).

Author Contributions

Conceptualization, C.W., P.B., and C.J.S.; Methodology, J.-S.B.; Investigation, C.W. and M.N.; Formal Analysis, R.S.-K.; Writing – Original Draft, C.W., P.B., and C.J.S.; Writing – Review & Editing, C.W., J.-S.B., M.N., R.S.-K., P.B., and C.J.S.; Funding Acquisition, P.B.; Supervision, P.B. and C.J.S.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and Algorithms
MATLAB 2016a	The MathWorks, Natick, MA, USA	https://uk.mathworks.com/products/new_products/release2016a.html
SPSS 20	IBM Corporation, Armonk, NY, USA	https://www.ibm.com/analytics/us/en/technology/spss/
FMRIB Software Library (FSL) 5.0	Oxford Centre for FMRI of the Brain (FMRIB), Wellcome Centre for Integrative Neuroimaging (WIN), Oxford, UK	https://fsl.fmrib.ox.ac.uk/fsl/fslwiki