Research ArticleNEUROSCIENCE

Dose-dependent effects of transcranial alternating current stimulation on spike timing in awake nonhuman primates

See allHide authors and affiliations

Science Advances  02 Sep 2020:
Vol. 6, no. 36, eaaz2747
DOI: 10.1126/sciadv.aaz2747


Weak extracellular electric fields can influence spike timing in neural networks. Approaches to noninvasively impose these fields on the brain have high therapeutic potential in neurology and psychiatry. Transcranial alternating current stimulation (TACS) is hypothesized to affect spike timing and cause neural entrainment. However, the conditions under which these effects occur in vivo are unknown. Here, we recorded single-unit activity in the neocortex in awake nonhuman primates during TACS and found dose-dependent neural entrainment to the stimulation waveform. Cluster analysis of changes in interspike intervals identified two main types of neural responses to TACS—increased burstiness and phase entrainment. Our results uncover key mechanisms of TACS and show that the stimulation affects spike timing in the awake primate brain at intensities feasible in humans. Thus, novel TACS protocols tailored to ongoing brain activity may be a tool to normalize spike timing in maladaptive brain networks and neurological disease.


It is increasingly recognized that local field potentials (LFPs) are not a mere by-product of synchronized neural activity but have an important function of controlling neural excitability (1, 2). This phenomenon, called ephaptic coupling, describes the coordinating influence of local electric fields or brain oscillations on neural spiking activity (3, 4). Spike timing is a key mechanism for neural encoding, communication in brain networks, and plastic changes (5, 6). Abnormal brain oscillations underlie a wide range of neurological and psychiatric disorders (7, 8), and tools to normalize pathological brain function have high potential for clinical applications.

It is possible to impose weak electric fields on the human brain in an intensity range comparable to endogenous fields using noninvasive technologies. One prominent method—transcranial alternating current stimulation (TACS)—applies weak oscillating electric currents through electrodes attached to the scalp (9). Cognitive neuroscientists increasingly use TACS in humans to manipulate attention (10), perception (11), and memory (12). Invasive intracranial measurements have demonstrated that intensities typically used in human studies, up to 2 mA (13), can achieve electric field strengths reaching up to 1 mV/mm (14, 15) in the human brain. A long history of in vitro studies and in vivo measurements in rodent models (1623) has established that these weak electric fields are sufficient to induce acute physiological effects. In vitro preparations (i.e., isolated cells or slices), however, lack crucial system-level properties. Rodent in vivo studies are limited due to the use of anesthesia, which affects calcium channels and neural dynamics. Furthermore, rodents have different brain anatomy and cytoarchitecture, which can lead to a different affinity to external TACS electric fields (24, 25). Thus, it is not clear how these findings translate to in vivo physiological effects in awake humans or nonhuman primates, where endogenous ongoing activity in complex oscillatory networks might either amplify or suppress the effects of weak externally applied electric fields.

Recent studies about the mechanisms of TACS are conflicting. One study using invasive electrophysiological measurements in surgical epilepsy patients did not find entrainment effects of TACS (26); however, only an indirect measure was used. Further, studies in anesthetized rats claimed that TACS intensities are too weak (27) or physiological responses arise from secondary stimulation effects of the peripheral nervous system (28). On the other hand, a recent study performing hippocampal recordings in awake nonhuman primates (29) showed that TACS can directly entrain single-neuron activity in a frequency-specific manner. However, TACS electric fields using the same stimulation montage are typically larger in nonhuman primates than in humans (14, 24) and findings in deep brain regions could be mediated by cortical effects experiencing higher electric fields. Thus, an investigation of direct stimulation effects in neocortex in a sub-to-milliamp intensity range is crucial to establish TACS as a meaningful tool to affect neural activity.

Here, we conduct dose-response measurements of single-unit activity in the pre- and postcentral gyri of awake nonhuman primates during TACS (Fig. 1). The stimulation is applied at 10 Hz, which is close to ongoing oscillatory activity measured with LFPs in this brain region [fig. S1, see also (30)]. Our approach overcomes several limitations of previous studies: First, we perform recordings in awake nonhuman primates. It is well known that anesthesia affects neural excitability, which can strongly influence the response to TACS. Further, a nonhuman primate model is more closely related to the human brain where gyrification is important for brain stimulation effects (31). Second, single-unit recordings are a direct electrophysiological measure of neural activity, which can be recovered in the presence of a stimulation artifact (23, 29). This is in contrast to previous studies measuring LFPs or scalp electroencephalogram recordings in humans, which are heavily affected by complex artifacts (26, 32). Third, simultaneous electric field recordings control the effective stimulation dose (33) and allow direct comparison to human applicable stimulation parameters.

Fig. 1 Experimental outline and data processing.

(A) A 96-channel microelectrode microdrive (maroon dots indicate different contacts at their original cortical location) was used to record from motor cortical regions (pre-motor and primary motor cortex). TACS was applied through two round stimulation electrodes (black) attached on the scalp over the left (L) and right (R) temples. Top to bottom is anterior to posterior direction. The brain orientation is identical for all panels. (B) Stimulation time course of different conditions. TACS was applied at 0.5, 1.0, or 1.5 mA for 2 min with resting periods before and after stimulation. (C) Recorded raw data of neural activity (orange) in the context of the TACS artifact and other noise sources (black). Signal filtering suppressed the TACS artifact, allowing spike extraction, sorting, and identification of time stamps for spike times with respect to the stimulation phase. (D) Overlaid spike waveforms before, during, and after TACS (from one exemplary neuron), demonstrating the consistency of spike waveforms in the presence of the stimulation artifacts.

Our findings demonstrate the effectiveness of TACS to affect neural spike timing and entrain spikes at the applied stimulation frequency in a dose-dependent fashion. With increasing electric field intensity, more neurons respond to stimulation. Last, by clustering interspike interval (ISI) changes across neurons, we identify two independent responses to TACS: increased burstiness and phase entrainment. The findings are a first mechanistic study of dose-response relationships of immediate effects of TACS in awake primates and provide important insight into the neurophysiological mechanisms of TACS.


Electric field recordings

Measured electric fields showed a dominant left-right orientation (Fig. 2A), which is expected given the location of the two stimulation electrodes. Median electric field strengths at the locations of the identified neurons were found to be 0.38, 0.77, and 1.15 mV/mm in subject A (Fig. 2A) and 0.43, 0.86, and 1.33 mV/mm in subject B for 0.5, 1.0, and 1.5 mA, respectively (Fig. 2B). These values are in good agreement with previous invasive measurements or modeling studies (14, 24). According to Ohm’s law, the electric field strength increases linearly with intensity, while the relative spatial distribution remains identical (fig. S2). Maximum electric field strength was found toward the lateral edge, which could be explained through the proximity of the stimulation electrode and gray matter–cerebrospinal fluid transition (34).

Fig. 2 Recorded electric field during TACS in the neocortex.

Electric field strength is color-coded, with red indicating maximum electric field strength. Right to left is anterior to posterior direction. The brain orientation is identical for all panels. (A) Electric field distributions in subject A on the brain surface for three different intensities (0.5, 1.0, and 1.5 mA). The distributions have the same orientation and relative spatial relationships and scale linearly with intensity. (B) Electric field distributions in subject B for three different intensities and for the shoulder control stimulation.

TACS dose-response relationship on spike entrainment

We evaluated the alignment of spike times with respect to the stimulation phase (“entrainment”). We found that several recorded neurons did spike closer in time to the TACS peak (0 degree) and less during the trough (Fig. 3). During both the TACS peak and trough, the electric field strength is maximal but with opposite direction. Depending on the orientation of the neuron with respect to the electric field, this is consistent with the idea of maximizing de- or hyperpolarization along the cortical neuron (35). With increasing stimulation intensity, spiking activity became increasingly phase-locked to the TACS peak times and suppressed during the TACS trough (Fig. 3A). This entrainment effect occurred immediately after stimulation onset and ceased after offset. Notably, our stimulation lasted only 2 min, which is likely not enough to induce lasting effects. While entrainment to TACS was pronounced in several neurons, it occurred only in a subset of all recorded neurons (Fig. 4). Thus, on the group-level, mean phase locking values (PLVs) increased only slightly during TACS (Fig. 4A). Identifying neurons based on their degree of phase entrainment (P < 0.01 in the Rayleigh test), we found that TACS at 0.5, 1.0, and 1.5 mA entrained 8.9, 17.6, and 26.5% of all neurons (n = 34), respectively. The responsive neurons were subjected to an electric field strength ranging from 0.2 to 1.2 mV/mm (see individual PLV-field data in Fig. 5). The generalized linear mixed effect model indicates a significant effect of stimulation intensity on the Rayleigh z value (F1,304 = 11.54, P = 0.0008) as well as on PLV (F1,304 = 7.13, P = 0.008). Together, 12 neurons (out of 34) demonstrated phase entrainment at some intensity level. Of them, nine cells showed a linear dose dependency: Two neurons were entrained starting from 0.5 mA, another two neurons were entrained starting from 1.0 mA, and five more neurons were entrained starting at 1.5 mA. Three remaining cells responded at the lower, but not higher, intensity. In the responsive neurons, the mean PLVs were increasingly enhanced with higher stimulation intensities (Fig. 4B; mean PLV = 0.11, 0.15, and 0.19, respectively). Out of 18 cases of phase locking, 11 show entrainment to the peak of the stimulation waveform (circular median = 0 ± 35 deg; see individual data in Fig. 6), 4 cases to the falling phase (90 ± 35 deg), and 3 more cases to the rising phase (270 ± 35 deg). These phase preferences were not dependent on the stimulation intensity (circular-linear correlation, P = 0.14). All neurons that responded at multiple intensity levels maintained consistent phase preference across intensities. During the pre- and post-stimulation periods, PLVs remained at the same nonsignificant level (mean PLV = 0.06, P > 0.01 in the Rayleigh test). In addition to commonly used PLV, we estimated a similar but unbiased statistical measure—pairwise phase consistency (PPC) (36). The findings are comparable and are shown for each neuron in table S1. The statistical model of PPC dependency from the stimulation intensity also captures the significant main effect (glmm F1,304 = 5.56, P = 0.019). Firing rates showed no consistent changes (Fig. 4, C and D; glmm F1,304 = 0.08, P = 0.78). We found no changes in the spiking activity in the shoulder control stimulation condition (fig. S3; all Rayleigh P > 0.1, mean PLV = 0.04).

Fig. 3 Dose dependency of single-unit entrainment.

(A) Time course of spikes (indicated as black dots for pre-/post-stimulation periods, and as orange dots during stimulation) over two stimulation cycles (x axis) over total recording time before, during, and after TACS (y axis). The data are shown for one exemplary neuron for three intensities (0.5, 1, and 1.5 mA). Pre-stimulation trials and post-stimulation trials (black) are arranged to a virtual stimulation signal with the same frequency (10 Hz) as during stimulation (orange). With stimulation onset, spike times cluster more toward the peak of the stimulation cycle (corresponding electric field shown in Fig. 2) compared to a more uniform distribution during rest. Below, peristimulus histograms (smoothed with a 10-ms Gaussian kernel) of spike times during stimulation (TACS waveform is shown with gray above the figure) indicate a nonuniform distribution of spike onsets. On the right side, polar plots of the phase during spike onset are shown. *P < 0.01, significant nonuniformity according to a Rayleigh test. (B) Same analysis as presented in (A) but for another exemplary neuron with higher spiking rate.

Fig. 4 Phase-lag value and spike rate changes across all neurons.

(A) Phase-lag values across all recorded neurons (n = 34) are shown for all conditions (before, during, and after stimulation) for three studied intensities (0.5, 1.0, and 1.5 mA). Individual dots indicate values for each neuron, and solid bars show the group mean. (B) The increase in PLV during stimulation is driven by a subset of neurons responding to TACS with increased entrainment. PLV values during stimulation are enhanced during all intensities and increasingly so for higher intensities. Bright blue and red highlight the cases from Fig. 3. (C) Spiking rate during TACS and rest across all conditions are shown for all recorded neurons and (D) for only responsive neurons. No consistent effects on spiking rate were observed during TACS. See fig. S3 for the same data during shoulder stimulation (control condition).

Fig. 5 Relationship between the electric field and PLV.

Linear relationship of the relative changes in PLVs, ΔPLV = (PLVstim − PLVpre)/PLVpre, to the local electric field strength (r.u., relative units). Each neuron (n = 34) is displayed during 0.5-, 1.0-, and 1.5-mA TACS. The responsive cases are highlighted in orange. Linear regression for the responsive cases: R2 = 0.39, F16 = 10.30, P = 0.005; for all cases: R2 = 0.01, F100 = 0.58, P = 0.45.

Fig. 6 Individual phase histograms of neural spikes during TACS.

Polar phase plots of the timing of neural spikes in individual “responsive” neurons relative to the stimulation waveform (10 Hz) during TACS for intensities of 0.5 mA (A), 1.0 mA (B), and 1.5 mA (C). All plots show a significant deviation from the uniform distribution according to a Rayleigh test (P < 0.01). Phase values are given relative to a cosine, i.e., 0° is the peak and 180° is the trough of the stimulation waveform. The number in the lower right corner of each polar plot indicates the unique neuron.

ISI changes during TACS

Investigating a larger range of possible neural effects of TACS apart from entrainment, we calculated ISI histograms during TACS and at rest. We determined changes in spike timing induced by TACS by calculating the difference between ISI distributions in the on- and off-stimulation conditions. We identified principled types of neural responses to TACS by calculating Newman’s modularity (37) based on the cross correlations of the ISI histograms. We identified two main classes of spike timing changes during TACS (Fig. 7 and fig. S4): increased burstiness (Fig. 7A, see also fig. S5) and phase entrainment (Fig. 7B, see also fig. S6). With increasing intensities, one class of neurons showed a shift toward shorter subsequent spike times during TACS compared to rest [n = 17 (9 + 8), 21 (10 + 11), and 18 (9 + 9) neurons (in subjects A + B) during 0.5, 1.0, and 1.5 mA, respectively]. This is not accompanied, however, by a significant increase in overall firing rate. A second class of neurons showed entrainment to TACS with increased subsequent spike times at 100 ms and an accompanying decrease in faster spikes [n = 17 (6 + 11), 13 (5 + 8), and 16 (6 + 10) neurons, respectively]. These two classes were found independently across all three stimulation intensities. The mean ISI histograms of the two classes had more extreme differences at every intensity level than one could expect from two randomly defined classes (cluster-based permutation test, P < 0.01). Notably, there is no statistical relationship between the ISI clustering and stimulation-induced changes in PLV (point-biserial rank correlation, P = 0.26). Analyzing the ISIs during shoulder stimulation, we found no meaningful clustering of neural responses (cluster-based permutation test, P > 0.1).

Fig. 7 Interspike intervals.

Two types of stimulation-induced ISI changes (TACS minus rest). Each subplot shows a mean ISI histogram of neurons that belong to a given class according to Newman’s modularity analysis. All mean histograms are statistically significant (cluster-based permutation test, P < 0.01). Differences in spike counts (n.u., normalized units) of ISI histograms across all neurons falling into the class are color-coded. (A) Class 1. An increase in burstiness (spikes following each at low interstimulus times) during TACS is observed. This effect increases for higher stimulation intensities. (B) Class 2. An increase in entrainment (enhanced ISIs at 100 ms, based on 10-Hz waveform) is visible during TACS. This is accompanied with a decrease in fast successive spikes. See figs. S5 and S6 for the ISI histograms of individual neurons during TACS and rest.


By conducting single-unit recordings in awake nonhuman primates, we showed that TACS affects spike timing in a dose-dependent fashion. With increasing electric field strength, more neurons entrained to the stimulation frequency. Further, we identified increased burstiness as a second type of neural response to TACS. These spike timing changes occurred in a sub-to-milliamp stimulation regime, which can be upscaled to feasible parameters in humans (24).

One question of key importance for translational efforts to use TACS to affect brain rhythms, e.g., to address oscillopathies in schizophrenia (7), is the intensity range at which physiological effects occur. If effective fields cannot be achieved in a safe and painless manner in humans, then the therapeutic potential of TACS is severely hampered. Our electric field recordings indicate that entrainment effects can occur at field strengths <0.5 mV/mm, which is achievable in humans for TACS intensities between 1 to 2 mA (14). However, physiological effects were more pronounced for higher intensities; thus, our results can serve as an important starting point for the development of principled TACS protocols in humans to achieve sufficient electric field strengths for robust physiological effects (38, 39).

Our data show that 10 to 20% of neurons respond to TACS with spike entrainment. Experiments in vitro (16, 35) and modeling studies (40) argue that the responsiveness of neurons depend on their morphology and orientation to the applied electric field as well as their dynamic biophysical properties (e.g., ion channel dynamics). One possibility is that the entrained neurons were optimally aligned with respect to the electric field. Another possible reason is that their cell type–specific biophysical properties rendered them more susceptible to stimulation. Further, the increase of entrained neurons for higher field strengths fits well with the observation of a linear membrane depolarization for weak electric fields (41). As a second neural response, we found increased burstiness of neurons. Burst firing of neurons has been observed in the presence of enhanced background activity in a low frequency range (42, 43). It is conceivable that, in our recording, TACS increased this background activity, thus resulting in the observed effect.

The time course of our recordings indicates that entrainment occurred and ceased immediately with stimulation on- and offset, respectively. This does not support theories suggesting a buildup phase or stimulation echo (44). However, our protocol used short stimulation durations (2 min) that were not intended to investigate the extent to which TACS elicits after-effects, such as due to spike-dependent plasticity (45, 46). Lasting effects of TACS have been reported in the literature following 10 or more minutes of stimulation (47, 48). Future studies systematically manipulating the TACS parameter range (frequency, electric field orientation, duration) will be important to identify the most effective stimulation parameters for immediate and lasting physiological effects. These studies will be facilitated by our increasing understanding of the underlying biophysics of TACS (14, 49) and sophisticated models to predict intracranial electric fields (34).

It is unlikely that the observed physiological responses were due to peripheral nerve stimulation or other peripherally mediated effects. The absence of entrainment during the control shoulder stimulation indicates a more direct effect of TACS on the brain. This is in line with another recent TACS study in awake nonhuman primates performing additional control measurements using local skin anesthesia (50). Notably, our measured strength of neural entrainment, PLVs, is noticeably higher than in the preceding study in anesthetized rodents that found no direct neural effect of stimulation (28). Thus, generally suppressed brain activity could complicate past findings. It is very likely that an awake brain is more susceptible to the modulation compared to an anesthetized state explored in the previous research (27, 28), thus explaining differences found in effective stimulation doses across studies.

Because TACS affects a large brain volume, it is likely that the stimulation also influenced lateral brain regions underneath the stimulation electrodes, possibly experiencing higher electric field strengths. Then, the lateral regions could have possibly propagated the effect into the recorded pre-central area. As it is not feasible to observe all neurons in the alive primate brain, the question of direct versus propagated neural effects of TACS remains a matter of active debate (51, 52). Still, one can infer practical conclusions from knowing the input to the brain (current intensity) and output (response of the observed neurons). Given that spike entrainment in this study occurred at different preferential phases, as was also observed previously (29), it is plausible that multiple mechanisms coexist. On another note, the electric field in our measurement brain area has a predominant left-right orientation along the neocortex that “connects” the two stimulation electrodes. Thus, the craniotomy above the recording area had little impact on the intracranial electric fields; otherwise, the field would have been oriented from the craniotomy inward to the brain, perpendicular to its surface.

In conclusion, our results show that weak extracellular electric fields can affect spike timing in vivo in awake nonhuman primates in a dose-dependent fashion. Enhancement or suppression of spiking activity during the peaks and troughs of stimulation suggest a neural membrane de- or hyperpolarization due to the applied weak extracellular fields. Entrainment effects were found in a subset of neurons that increased with stimulation intensity. Another robust neural response—increased burstiness—has not been considered so far in the TACS literature. This suggests the need for further systematic studies to examine TACS mechanisms in vivo. The found physiological responses occurred in a sub-to-milliamp intensity range potentially transferable to humans, thus providing an electrophysiological grounding for efforts to interact with human brain oscillations in a noninvasive and safe manner.


Surgical procedures

All procedures were approved by the Institutional Animal Care and Use Committee of the University of Minnesota and complied with the U.S. Public Health Service policy on the humane care and use of laboratory animals. Two adult female nonhuman primates (Macaca mulatta) were used in this study (subject A: female, 10 kg, 18 years; subject B: female, 12 kg, 21 years). The animals were implanted with a head post and cephalic chamber positioned midline over motor cortical areas. In subject A, all implant materials were made of magnetic resonance imaging (MRI)–compatible, nonconductive polyetheretherketone (PEEK) and ceramic except for a single titanium bone screw; in subject B, the implant materials including headpost, chambers, and bone screws were made of titanium. After surgery in both animals, a 96-channel microdrive with independently moveable tungsten microelectrodes was fixed to the cephalic chamber (Gray Matter Research). Confirmation of the chamber placement and electrode positions was done by co-registering preoperative 7-Tesla T1 MRI and postoperative computed tomography (CT) scans using the Cicerone software package (53).

Electrophysiological recordings and transcranial current stimulation

Neurophysiological data were collected using a TDT workstation (Tucker-Davis Technologies) operating at ∼25-kHz sampling rate. A large dynamic range (±500 mV) and 28-bit resolution allowed recording of electrophysiological data during TACS, fully capturing both neural data and stimulation artifacts without data loss (e.g., due to amplifier saturation). In subject A, the common reference for the electrophysiological recordings was a titanium bone screw positioned anterior to the recording chamber; in subject B, the common reference was the cephalic chamber. All data were collected during an awake, resting state while the animal was seated in a primate chair that kept the head facing forward. TACS was applied through two round stimulation electrodes [3.14 cm2, Ag/AgCl with conductive gel (SignaGel)] attached over the skin of the left and right temple (Fig. 1A). This electrode montage was chosen to minimize effects of the implant and skull defects on the TACS electric fields (54) by applying currents in an orthogonal direction to the recordings (15). In addition, performing concurrent electric field recordings (see below) allowed us to directly measure and control the biophysical stimulation condition. This ensured that results would be comparable to stimulation parameters commonly used in human TACS studies.

In total, we had six main experimental blocks (three intensities in two subjects). At the beginning of the experiments, we recorded a rest condition without stimulation for several minutes. Then, TACS was applied at a frequency of 10 Hz for 2 min with 5-s ramp up and down of the current amplitudes. Stimulation intensities applied were 0.5, 1.0, and 1.5 mA (peak-to-baseline, Fig. 1B). Before and after each TACS run, resting condition was recorded for several minutes. All neurons were recorded at the same time. In addition, a peripheral control stimulation block was conducted where TACS was applied at 10 Hz and 1 mA intensity through two electrodes attached over the shaved right upper arm in subject B.

Data analysis

One key issue to evaluate neurophysiological TACS effects measured with electrophysiology is the recovery of neural activity in the presence of a stimulation artifact. TACS induces a narrow band sinusoidal signal in the range of several millivolts (14, 55) at the applied stimulation frequency. Because of the presence of harmonics and interactions with other physiological signals, such as heartbeat and breathing (32, 5658), a full recovery of LFP data in the spectral range of interest is challenging. Spikes are in a frequency range much higher than the TACS signal; thus, filtering can suppress most of the stimulation artifact (see Fig. 1C for graphical outline of data processing). Even more important, spikes have a characteristic waveform that can identify them in the presence of other superimposed signals (Fig. 1D). Because of these considerations, we focused our analysis on spiking activity during TACS.

Spike sorting

Within the 96-channel microdrives, neuronal spike recordings were collected from 15 (subject A) and 11 (subject B) microelectrodes that were inserted in pre- and postcentral gyrus regions. These wide-band recordings were then analyzed offline using custom software developed in MATLAB (MathWorks) and Offline Sorter (Plexon). Raw signals were bandpass-filtered between 0.3 and 10 kHz in subject A and 0.3 and 5 kHz in subject B, and single units were isolated and sorted using principal component and template-based methods in Offline Sorter. In subject B, TACS artifacts were larger than in subject A (likely due to the difference in recording ground). As such, a comb filter was also used to remove artifact harmonics that bled into the neural spike recording spectrum (1 to 5 kHz). We identified 38 neurons (n = 17 for animal A and n = 21 for animal B), of which four where excluded due to very sparse firing (<0.25 Hz) not resulting in a sufficient number of spikes in the studied time period (final n = 15 for animal A and n = 19 for animal B).

Electric field recording

All 96 contacts within the microdrive were used for estimating the electric fields during TACS. In both subjects, the electrodes covered a space of approximately 16 mm × 16 mm on the surface plane and approximately 12 mm in depth, forming a three-dimensional (3D) cloud with millimeter-to-submillimeter spacing between the closest contacts. This enabled us to accurately capture all directions of the electric field. To compute the electric field, the recorded data were first bandpass-filtered around the applied stimulation frequency (10 ± 1 Hz) using a forward-reverse, zero-phase, second-order Butterworth IIR filter to derive the TACS voltage signal (14). This resulted in a 3D voltage distribution at the location of the recording electrodes. Electric field E in the recording area was then estimated by linearly interpolating the recorded voltages in a regular 3D grid and calculating the numerical derivative across the grid points as i(∂V/∂x) + j(∂V/∂y) + k(∂V/∂z), where V is the voltage potential, and i, j, k are the unit vectors in x, y, z direction (Fig. 2).

Dose-response of spike entrainment to stimulation phase

In a first analysis, we evaluated dose-response effects of spike entrainment to TACS. For this, we determined the relative onset of spike times with respect to the TACS signal and computed peristimulus histograms. These show the number of spikes per time period of the stimulation. Deviation from uniformity indicates a preference for spiking relative to distinct stimulation phases. To compare the stimulation-induced changes to the resting “control” condition, we determined spike times at rest (before and after stimulation) relative to a virtual sinusoid with the same frequency as during TACS. In addition, we computed phase histograms indicating the preferred TACS phase during spiking (0 degree equals peak of stimulation). Further, we computed firing rates and PLVs before, after, and during stimulation for all neurons. PLV = abs(mean(exp(iθ))), where θ is the phase of the stimulation waveform at which the spike occurs. We identified neurons as responding to TACS based on a Rayleigh uniformity test (Rayleigh z value > 4.6, P < 0.01) using the circular statistics toolbox (59). These analyses were performed and compared across all three stimulation intensities (0.5, 1.0, and 1.5 mA) and during the shoulder stimulation (1.0 mA).

To test the statistical significance of changes in Rayleigh z values, PLVs, and firing rates, we implemented three generalized linear mixed effect models (glmm). Each outcome variable had 306 observations (34 neurons × 3 intensities × 3 conditions, which are the periods before, during, and after stimulation). In each model, the stimulation intensity (0.5, 1.0, or 1.5 mA) was the fixed effect factor, and the presence of stimulation (on or off), the subject ID (subject A or B), and the membership of the neuron to an ISI cluster (first or second cluster, see the next section) were the random effect factors to account for possible interindividual variability in the neural responses. The model parameters were estimated using a restricted maximum pseudo-likelihood procedure with a log link function. All F values, degrees of freedom, and P values are reported in the main text.

ISI changes during TACS

While entrainment is one of the key postulated mechanisms for TACS (39), it is conceivable that other changes in spike timing can occur due to stimulation. To study this, we computed ISI histograms during rest and during TACS. ISIs indicate the time difference between two subsequent spikes (e.g., 0.1 s if spike n is followed by a spike n + 1 with a 100-ms delay). ISI histograms plot the ISI of spike n against the ISI of spike n + 1. ISI histograms can give insights into the temporal spiking patterns of neurons (e.g., indicate burstiness if spikes are followed by fast subsequent ones or spiking periodicity when spikes occur at similar time delays). We computed ISI histograms during TACS and resting state conditions and estimated their difference (TACS – rest). To identify independent classes of neural responses, we computed a 2D cross-correlation matrix from the ISI histograms, which can be treated as a weighted graph. Using Newman’s modularity (37) implemented in the Brain Connectivity Toolbox (60), we identified two main classes of TACS responses. Newman’s modularity is a community detection algorithm; it is searching for clusters in the data that have rich connections (in our case, high cross-correlation) within the cluster and poor connections (low cross-correlation) outside the cluster. All ISI difference histograms were then averaged within each class and intensity to create a mean class histogram. The statistical significance of the mean class histograms was assessed using a cluster-based permutation test (12,500 permutations). For this, we shuffled the class labels of the individual histograms, computed dummy mean class histograms (for each stimulation intensity), and compared the critical values of the real mean class histograms with those of the dummy mean class histograms. We defined the critical value as the maximum number of adjacent nonzero elements in the histogram after binarizing with a 90% intensity threshold.


Supplementary material for this article is available at

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.


Acknowledgments: We thank K. Wilmerding and S. Nebeck for assistance with the experimental preparations. Funding: Research reported in this publication was supported by the University of Minnesota’s MnDRIVE Initiative, in part by the NIH (P50-NS098573, R01-NS058945, R01-NS094206), a predoctoral fellowship (NINDS, F31-NS108625), and Engdahl Philanthropic Donation. Author contributions: L.J., I.A., and A.O. designed the experiments; L.J., J.K., A.D., and A.O. collected the data; L.J. and M.J. supervised the data collection; L.J., I.A., J.K., A.D., Y.Y., and A.O. processed MR and CT data and identified electrode locations; L.J., I.A., M.J., and A.O. analyzed the electrophysiological data; M.J. and A.O. supervised the data analysis; L.J. and I.A. prepared the figures; L.J., I.A., M.J., and A.O. interpreted the results; A.O. wrote the paper with critical input from the coauthors; all authors reviewed the paper. Competing interests: The authors declare that they have no 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. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article