Abstract
Transcranial direct current stimulation (tDCS) has emerged as a potentially safe and effective brain stimulation modality that alters cortical excitability by passing a small, constant electric current through the scalp. tDCS creates an electric field that weakly modulates the membrane voltage of a large number of cortical neurons. Recent human studies have suggested that sine-wave stimulation waveforms [transcranial alternating current stimulation (tACS)] represent a more targeted stimulation paradigm for the enhancement of cortical oscillations. Yet, the underlying mechanisms of how periodic, weak global perturbations alter the spatiotemporal dynamics of large-scale cortical network dynamics remain a matter of debate. Here, we simulated large-scale networks of spiking neuron models to address this question in endogenously rhythmic networks. We identified distinct roles of the depolarizing and hyperpolarizing phases of tACS in entrainment, which entailed moving network activity toward and away from a strong nonlinearity provided by the local excitatory coupling of pyramidal cells. Together, these mechanisms gave rise to resonance dynamics characterized by an Arnold tongue centered on the resonance frequency of the network. We then performed multichannel extracellular recordings of multiunit firing activity during tACS in anesthetized ferrets (Mustela putoris furo), a model species with a gyrencephalic brain, to verify that weak global perturbations can selectively enhance oscillations at the applied stimulation frequency. Together, these results provide a detailed mechanistic understanding of tACS at the level of large-scale network dynamics and support the future design of activity-dependent feedback tACS paradigms that dynamically tailor stimulation frequency to the spectral peak of ongoing brain activity.
Introduction
Successful brain stimulation paradigms either disrupt pathological activity or enhance physiological processes. Clinically, there is a demand for advanced neuromodulation paradigms based on growing evidence that neuropsychiatric illnesses such as schizophrenia are “network disorders.” Such conditions are fundamentally mediated by aberrations in the synchronization and organization of oscillatory cortical network activity (Uhlhaas and Singer, 2012) and are likely mediated by deficits in network connectivity (Zalesky et al., 2011; Samartzis et al., 2013). The development of stimulation paradigms capable of modulating specific frequency bands has been hampered by the lack of mechanistic understanding of how neuronal networks respond to applied stimulation.
This is particularly true for transcranial current stimulation (TCS), an old, noninvasive stimulation modality that is undergoing a renaissance due to (1) small, yet promising early trials for cognitive enhancement and treatment of neurological and psychiatric conditions (Paulus, 2011; Brunoni et al., 2012) and (2) a growing understanding of how weak electric fields interact with networks of neurons (Deans et al., 2007; Radman et al., 2007; Fröhlich and McCormick, 2010; Ozen et al., 2010; Reato et al., 2010). Constant waveforms [transcranial direct current stimulation (tDCS)] either suppress or enhance activity as a function of electrode polarity (Nitsche and Paulus, 2000). Recent studies have demonstrated that transcranial alternating current stimulation (tACS) can induce frequency-specific effects on brain dynamics measured by electroencephalography (EEG; Kirov et al., 2009; Kanai et al., 2010; Moliadze et al., 2010; Zaehle et al., 2010; Zaghi et al., 2010) and behavioral tasks. Specifically, this targeted stimulation paradigm has been applied to modulate vision (Kanai et al., 2008; Laczó et al., 2012), motor function (Antal et al., 2008; Pogosyan et al., 2009; Feurra et al., 2011b; Brignani et al., 2013), somatosensation (Feurra et al., 2011a; Wach et al., 2013), and higher-order cognitive function (Polanía et al., 2012; Sela et al., 2012; Brittain et al., 2013). The effects of weak global perturbations on larger cortical networks, such as the estimated millions of neurons that are affected in human TCS, have remained unstudied. Filling this gap in knowledge is important because small network models often exhibit fundamentally different dynamics than their larger counterparts.
Here, we use large-scale simulations of cortical networks to determine whether and how weak global perturbations modify the macroscopic behavior of large-scale networks. We hypothesized that (1) simulated tACS matched in temporal structure to spontaneous network activity will stabilize and enhance spatiotemporal activity motifs that occur in the absence of stimulation, and (2) spatial patterning of activity in these large networks plays an important role in mediating TCS effects. We verified the modulatory effects of simulated TCS by performing in vivo electrophysiological recordings during TCS in anesthetized ferrets. This model species is ideal for studying TCS because animals exhibit gyrencephalic brains similar to humans, and therefore the effect of TCS is nonuniform across neurons. Gaining a mechanistic understanding of how weak global current perturbations interact with large-scale network dynamics will allow for the design of stimulation waveforms that may be safer and more efficient.
Materials and Methods
Intrinsic dynamics of model neurons.
The computational model used in this study consists of Izhikevich model (Izhikevich, 2003, 2004; Izhikevich and Edelman, 2008) pyramidal cells (PYs) and inhibitory interneurons (INs). As any other mathematical model of biological processes, these neuronal models are abstractions. This model was chosen for this study due to its realistic dynamics of spike initiation despite the low number of free parameters that guide the model neuron dynamics (Izhikevich, 2007). The intrinsic firing behavior of the model cell types used here was previously demonstrated to be virtually indistinguishable from neuronal recordings in vitro (Izhikevich, 2003). Briefly, the neuronal dynamics are described by two coupled differential equations with four parameters (a–d) that determine intrinsic excitability. The update rule for membrane potential, V′, of each neuron is as follows:
where V is the previous membrane voltage, Δt = 0.1 ms is the time step, EAMPA = 0 mV is the AMPAergic excitatory reversal potential, EGABA = −80 mV is the GABAA inhibitory reversal potential, GEX and GIN are the sum of all afferent excitatory (gPY) and inhibitory (gIN) conductances, ITCS and INoise are currents injected to model TCS and to induce spontaneous background noise, and u is the slow-recovery variable. For PYs, the values of parameters a (time scale of recovery) and b (sensitivity of recovery) were set to 0.02 and 0.2, respectively. To model regular spiking, intrinsically bursting and chattering cells in the PY layer, the reset potential parameter, c, ranged from −65 to −50 mV, and recovery after spike, d, from 6 to 8 across PYs. Both of these values are drawn from generalized Pareto distributions (μ = −50, σ = −30, ξ = −2, median = −61.26 mV for parameter c; μ = 6, σ = 4, ξ = −2, median = 7.50 for parameter d) that bias the parameter values such that regular spiking cells are the most common cell type in the PY layer. For the INs, the parameters c and d are fixed to −65 mV and 2, respectively. To model both fast-spiking and low-threshold neurons, parameters a and b are drawn from uniform distributions (0.02–0.1 and 0.2–0.25, respectively) such that either cell type is equally likely to occur in the IN layer.
Model of synaptic dynamics.
Synapses are modeled by steps in synaptic conductance in the case of a presynaptic action potential, and by exponential decay otherwise. All synapses of a given type are grouped per neuron for computational efficiency. A similar approach has been successfully used before (Fröhlich et al., 2008). The respective conductances are updated as follows:
where GEX and GIN are the respective total conductances at the occurrence of the last presynaptic action potential from the corresponding cell type, τEX = 2 ms and τIN = 3 ms are the decay time constants, and Δtpsp is the time since the last presynaptic action potential. PY–PY connections exhibit short-term synaptic depression (Tsodyks and Markram, 1997) with a single depression variable D (D = 1 indicates no depression; D = 0 indicates synapse fully depressed) that exhibits a single-exponential recovery time course (τD = 300 ms). PY–PY synaptic gPY–PY strength is given by the following:
where GPY–PY denotes the synaptic strength in the absence of depression (D = 1) and D is updated at each presynaptic action potential for all PY–PY synapses formed by that neuron:
where r = 0.6 represents the depression coefficient, indicating the fraction of synaptic resources available immediately after vesicle release caused by an action potential.
Network topology.
A 200,000 neuron network model was created by arranging both the PY and IN layers into two-dimensional (2-D) squares with 160,000 (400 × 400) and 40,000 (200 × 200) neurons, respectively. The choice of this network size was motivated by a tradeoff between computational cost and biological plausibility. A large, two-dimensional architecture has previously been shown to give rise to more sophisticated network dynamics than smaller, one-dimensional structures (Bazhenov et al., 2008). The topology of the local excitatory connections exhibited a sparse local connectivity scheme with each PY connecting to a random 40% subset of 121 cells in its surrounding 11 × 11 grid of PYs (GPY−PY = 0.06). Synaptic inhibition exhibited global random connectivity both for PY–IN excitation (GPY–IN = 0.0001; each PY connected to 25 random INs) and feedback IN–PY inhibition (GPY–IN = 0.0002; each IN connected to 49 random PYs). The choice of synaptic connectivity topology was motivated by the goal of reducing the complexity of network dynamics to the point where synaptic inhibition contributed little to the spatial dynamics but rather acted as a global negative-feedback mechanism, counteracting the positive feedback of the recurrent excitatory connectivity. This way, without extensive parameter tuning, slow rhythmic activity occurred in the network through the temporal interaction of this positive- and negative-feedback mechanism. The strong local connectivity in the PY layer allowed for more realistic spatiotemporal dynamics that lacked pronounced large-scale synchrony typically seen in smaller network models that oscillate.
Noise and stimulation currents.
To break the symmetry of the model and to create more realistic temporal variability, all cells received a current injection (INoise) that was the sum of (1) a constant current injection ranging from 0 to 1.5 drawn from a generalized Pareto distribution (μ = 1, σ = −3, ξ = −3, median = 0.1895) to create a spontaneously active subset of PYs and (2) a variable current with a random value drawn at every time step (uniform distribution from 0 to 2 for PYs, and 0 to 1.5 for INs). TCS was modeled with a small current injection (ITCS; range 0.00–0.13; corresponding to 0–13 pA, resulting in an average membrane voltage depolarization of 0–130 μV) into PYs, which are more susceptible to applied electric fields than inhibitory interneurons because of PY's elongated somato–dendritic axis (Tranchina and Nicholson, 1986; Lopez et al., 1991; Radman et al., 2009).
Data analysis and presentation.
Frequency analysis was performed by convolving the excitatory network signal (percentage of PYs firing as a function of time) with a family of Morlet wavelets of 1–10 Hz (Goupillaud et al., 1984). To quantify the spatial activity of the network, an algorithm was developed in MATLAB (Mathworks) to count the number of local spots of PY activity (“hotspots”). Two-dimensional activity maps were constructed that color code the average firing rate (over a 10 ms bin). Using the MATLAB Image Processing Toolbox (Mathworks), each resulting frame was converted to a binary image to separate active from quiet neurons. To account for noise in the resulting image, the bwareaopen function was used to remove connected neighborhoods of ≤10 activated neurons. The remaining eight connected neighborhoods of activated neurons were identified and counted using the bwconncomp MATLAB function (see Fig. 4A, activity hotspots identified by algorithm, each shown in a unique color). Spectral power is shown on a normalized scale from 0 to 1 and is referred to as relative power.
Ferret in vivo electrophysiology.
To determine the effect of tACS weak electric fields on multiunit action potential firing in vivo, extracellular electrophysiological recordings were performed in two adolescent ferrets (Mustela putoris furo, female, 15–17 weeks old). As these animals had not yet reached sexual maturity, possible alterations in regional connectivity or excitability induced by estrous cycle phase can be excluded. Since the effect of TCS depends on the relative angle between the somato–dendritic axis of the neurons and the direction of the applied electric field (Tranchina and Nicholson, 1986), we performed these experiments to establish that the perturbation dynamics studied in our models also occur in gyrencephalic brains. Briefly, animals were deeply anesthetized with an initial intramuscular injection of ketamine (30 mg/kg) and xylazine (1–2 mg/kg), intubated (to allow delivery of isoflurane anesthetic in 100% medical oxygen), and artificially ventilated (10–15 cc/stroke, 40–50 breaths/min). The deep anesthetic plane was maintained throughout the experiment with a combination of 0.5–1.0% isoflurane and xylazine (0.015 ml/h) administered intravenously with 5% dextrose lactated Ringer's solution (4.25 ml/h). At the time of stimulation and recording, ketamine from the initial anesthesia induction was no longer active because multiple hours had elapsed from the time of anesthesia induction to recording (Quibell et al., 2011). Animals were positioned on a heating blanket, and their physiology (electrocardiogram, pulse oxygen level, end-tidal CO2, and rectal body temperature) was continuously monitored to ensure continued deep anesthesia. End-tidal CO2 was between 30 and 50 mmHg for all experiments (Kohn, 1997). Surgical procedures consisted of an initial midline incision of the scalp, removal of underlying soft tissue, and a small (<2 mm) craniotomy over the lateral gyrus (primary visual cortex, 3 mm anterior and 9 mm lateral to bregma). After durectomy, the brain was covered with 4% sterile agar in physiological saline and a 16-channel linear recording depth probe (100 μm contact site spacing; NeuroNexus) was slowly advanced into cortex with a micromanipulator (Narishige) to record simultaneously from all cortical layers. For the electric reference of the neuronal recordings, a silver chloride wire was placed between anterior skull and retracted soft tissue, and was held in place with agar. The correct depth of the linear recording probe was determined by small deflections of the local field potential (LFP) at superficial electrode recording sites and larger deflections of the LFP at deeper electrode recording sites. Unfiltered signals were first amplified with MPA8I head stages with gain 10 (Multichannel Systems), then further amplified with gain 500 (model 3500, A-M Systems). Signals were high-pass filtered (300 Hz cutoff) to extract multiunit firing without contamination from the stimulation, digitized at 20 kHz (Power 1401, Cambridge Electronic Design), and digitally stored using Spike2 software (Cambridge Electronic Design). TCS was provided through two bilaterally positioned silver chloride electrodes (1 × 2 cm) that spanned the exposed anterior–posterior axes. Current was supplied by a voltage-controlled current source (peak amplitude, 2 mA; stimulation frequencies, 0.5 to 2.5 Hz in steps of 0.5 Hz; A395 Linear Stimulus Isolator, World Precision Instruments). Resistance across stimulation electrodes was constant throughout the experiment. Stimulation epochs lasted 30 s, were interleaved with 60 s of control before and after stimulation, and were randomized in counterbalanced order for stimulation frequency (fs) applied. Baseline spontaneous activity was constant before stimulation and recovered completely after stimulation. All recordings were conducted in a dark and quiet room.
Results are based on 480 trials for each stimulation frequency across two animals. MU activity was extracted by applying a negative threshold of −3 × SD of the recorded signal, and the spectrograms were determined by wavelet transformation (same as described above for model data) applied to the spike train smoothed with a Gaussian kernel (SD, 50 ms). Changes in spectrum were determined by averaging the spectrograms over time and subtracting the spectrum during control from the spectrum during TCS. Enhancement was then determined as the increase in power of the multiunit activity at the stimulation frequency. Significance was determined by Wilcoxon rank-sum test (significant for p < 0.05). The top 50% of channels on the depth probe were assigned to be “superficial,” and the bottom 50% of channels on the depth probe were assigned to be “deep.” Phase histogram of activity entrainment was determined by binning spike times by phase of the applied stimulation. Entropy was used to compute the nonuniformity of the distributions (i.e., how large the entrainment effect was). All animal experiments were conducted under an institutionally approved animal protocol and exceeded the relevant guidelines set forth by the NIH and the U.S. Department of Agriculture.
Results
Model of large-scale spatiotemporal oscillations in cortex
We built a large-scale computational model of a cortical network to understand how weak global perturbations by TCS alter macroscopic network dynamics. Specifically, we studied the effect of tACS on endogenously oscillating cortical networks. Our model consisted of two layers of regular spiking excitatory PYs and fast-spiking INs, respectively (Fig. 1A, top four traces, sample voltage traces for PYs, bottom trace, IN). Cortical network dynamics are shaped in part by the relative amount and timing of synaptic excitation and inhibition (Haider and McCormick, 2009). In our model, we chose the strength of recurrent synaptic excitation and feedback inhibition such that spontaneous activity in the network of PYs triggered network-wide transient epochs of activity (UP states) that were terminated by synaptic inhibition provided by the layer of INs. UP states were followed by epochs of quiescence during which the excitatory synapses recovered from synaptic depression, which in turn enabled the generation of the next UP state (Fig. 1B, color map of membrane voltages, C, percentage of activation of PYs and INs as a function of time). In essence, the relatively loose coupling between excitation and inhibition and the local connectivity resulted in rhythmic (∼3 Hz), albeit somewhat irregular, overall macroscopic activation of the network (Fig. 1D, relationship between activity of PY and IN layers). Due to the local connectivity scheme of the recurrent excitation, the network exhibited spatiotemporal activity patterns characterized by initiation of activity as a local hotspot that then expanded and propagated outward (Fig. 1E, two-dimensional map of PY layer membrane voltage during an UP state). Given these rhythmic spatiotemporal activity dynamics, our model allowed us to (1) study the mechanisms of tACS both in space and time, and (2) use parametrically varied stimulation waveforms to identify optimal stimulation paradigms based on the resulting mechanistic understanding of tACS at the network level (rational design).
Global oscillatory activity in computational model of cortical network. A, Sample membrane voltage traces of individual PYs and INs. B, Activity map of the PY (top) and IN (bottom) layer of the network shows periods of activity (UP states) and quiescence (DOWN states; cool colors, hyperpolarization; warm colors, depolarization). Two-dimensional layers of neurons flattened to one-dimensional representation. C, Percentage of neurons firing in the PY (top) and IN (bottom) layer of the network (4 s duration). D, Phase-plane representation of PY and IN activity demonstrates that a structured, temporal relationship of the activity between the excitatory and inhibitory networks gave rise to a global oscillation pattern. E, Two-dimensional maps of PY activity during a network UP state (membrane voltage color coded as above). Network activity originated in hotspots and subsequently expanded through lateral excitatory PY–PY connections.
tACS outperformed tDCS
Today, most TCS studies use direct current waveforms. Given the ubiquity of oscillations in cortex, tACS has recently emerged as a potentially beneficial alternative; we compared the effects of these two stimulation paradigms on the spatiotemporal dynamics of our model cortical network. The two stimulation paradigms were matched in amplitude (equivalent current set to 9 pA), and the tACS frequency was matched to the intrinsic frequency of the network of 3 Hz (Fig. 2A, left, network activity of both layers of networks stimulated by tDCS, and right, stimulated by tACS). We found that tDCS had only a very limited effect on the overall oscillatory network dynamics, whereas tACS was very effective at entraining the network at the frequency of the applied stimulation (Fig. 2B, left, tDCS, average relative power = 0.29 at 3 Hz, right, tACS, average relative power = 0.99 at 3 Hz).
Comparison of direct (tDCS, left) and alternating (tACS, right) current stimulation. A, Overall network activity of PY and IN layers with applied stimulation waveforms (left: tDCS; right: tACS). An fs = 3 Hz was chosen to approximately match the intrinsic network oscillation frequency. Stimulation (9 pA, modeling TCS) was applied to both networks. B, Spectrograms of the two simulations demonstrate the dramatic power increase at the simulated frequency during stimulation for the sinusoidal stimulation (right; average relative power at 3 Hz was 0.99 at stimulation onset) compared with the direct current stimulation (left; average relative power at 3 Hz was 0.29 at stimulation onset).
We first established whether this pronounced enhancement depended on the onset phase of the sine-wave stimulation relative to the endogenous network oscillation. Specifically, we conducted a series of simulations applying tACS at 3 Hz and amplitude of 9 pA, and observed changes in the network as the onset phase (defined as the starting phase of the applied stimulation) was varied in 20 increments from 0 to 2 π. We found that, independent of the onset phase, the stimulated network was maximally entrained after a few seconds. However, the time it took to entrain (entrainment time was defined as the time needed to reach 90% of peak steady-state power at the stimulation frequency) varied with onset phase (Fig. 3A; relative power of PY activity at the stimulated frequency of 3 Hz; warmer line colors indicate onset phases closer to π). From an onset phase of 0 to just lower than π, the entrainment time increased near-monotonically from 1.6 to 2.1 s (Fig. 3B, entrainment time as a function of onset phase). This range of onset phases corresponded to a depolarizing stimulus being applied at the onset of stimulation when the network was near the end of an UP state (Fig. 3B, inset 1, PY activity for a network which received stimulation with an onset phase of 2.83). The depolarizing stimulation triggered a new, yet lower-amplitude UP state out of phase with the network activity, which disrupted the intrinsic network oscillation and resulted in longer entrainment times. For starting phases closer to π (i.e., that provided shorter duration of depolarization), the new out-of-phase UP states occurred with even lower amplitudes. As a result, these networks required more cycles to lock to the applied stimulation.
Time to phase lock after the onset of tACS. A, Change in power of networks stimulated at 3 Hz with 9 pA starting at different onset phases. The line color indicates the onset phase of the stimulation waveform at stimulation onset for that trial (increasing onset phase with warmer colors). All onset phases eventually entrained the network. B, Entrainment time (defined as 90% of peak steady-state power at the stimulation frequency) as a function of the onset phase of stimulation. The entrainment time was minimal for the network stimulated with the onset phase of stimulation at π. Insets 1 and 2 show the PY network activity for networks stimulated with the onset phase of stimulation set to 2.83 and 3.14 radians, respectively.
However, once the starting phase of stimulation was set to π, the entrainment time was minimal at 0.6 s (Fig. 3B, inset 2, PY activity for a network that received stimulation with an onset phase of π). In this case, the hyperpolarization provided at stimulation onset co-occurred and accelerated the transition into the next DOWN state. Furthermore, for onset phases of stimulation greater than π (i.e., that provided shorter duration hyperpolarization), a depolarizing stimulus triggering an UP state occurred sooner than it would have without stimulation. This disrupted the phase of network oscillations, and, as a result, the time to entrain again increased nearly monotonically for stimuli with starting phases larger than π. Together, these results demonstrate that the choice of onset phase of tACS determined the time for the network to entrain but had little effect on the steady-state entrainment.
Increase in hotspots as the network became more regular
To better understand the mechanisms of entrainment, we studied the effects of tDCS and tACS on the generation of activity hotspots that mediated the onset of UP states. We quantified the number of hotspots present in the PY layer as a function of time (Fig. 4A, during control periods, PY hotspots are color coded). During control periods (no stimulation), hotspots were nearly always (98% of the time) present. With tDCS, the fraction of time with hotspots was reduced to 78% (Fig. 4B, PY hotspot and network activity for PYs of both tDCS, top two plots, and tACS, bottom two plots, stimulated networks). However, with tACS, bimodal peaks occurred in the time course of hotspots during UP states interleaved by pronounced DOWN states (hotspots present for only 48% of the time).
The number of activity sites in the network increased as the network became more synchronized. A, Sites of activity (hotspots) were identified by counting groups of neurons with high firing rates. B, Top two plots, Network activity of the PY layer of a network stimulated using 9 pA tDCS and the number of activity sites. Bottom two plots, Hotspot activity for a network stimulated using tACS at 3 Hz and 9 pA. Unstimulated and tDC-stimulated conditions resulted in uniform hotspot activity, but tACS caused formation of bimodal peaks in the time course of hotspots. C, Plots of the relationship between network activity and the number of hotspots for unstimulated (left), tDCS (middle), and tACS (right) networks. The network was only fully activated with tACS; the maximum number of hotspots was seen just before or after an UP state.
UP states during control periods were typically triggered by a single hotspot that expanded outward by local excitatory connections. However, this expansion of the hotspot was slower than the collapse of its core. As a result, the maximum percentage of PYs active was <85% (Fig. 4C, number of hotspots as a function of percentage PY activity, left, control, middle, tDCS). In the case of tACS (Fig. 4C, right), UP states were initialized by multiple hotspots appearing nearly simultaneously. Hotspots appearing during entrained UP states coalesced quickly, resulting in 98% of PYs being active during the peak of an UP state (single hotspot that covers entire PY network). Toward the end of an UP state, PY activity collapsed and fragmented into many individual hotspots in the PY network.
These results describe how tACS entrained the network to the stimulation frequency by modulating the spatiotemporal structure of the network. When stimulated with tACS, a larger number of hotspots was generated by the network at the onset of an UP state compared with unstimulated or tDCS-stimulated networks. This high density of hotspots converged more quickly, activating the entire network so that all neurons fired during an UP state.
Recovery of synaptic depression played an important role in entrainment by tACS
To understand the role of negative and positive modulation of the membrane voltage with tACS, we applied decomposed stimulation waveforms that had either only the positive half-wave (i.e., “depolarizing-only”) or the negative half-wave (i.e., “hyperpolarizing-only”), with the amplitude of the other half-cycle set to zero. We found that hyperpolarizing-only stimulation entrained the network just as effectively as full-wave tACS. This was not the case for the depolarizing-only stimulation, which exhibited a pronounced, yet reduced, effect on the ongoing oscillation dynamics compared with full-wave tACS (Fig. 5A, PY network activity for depolarizing-only stimulation, left, with average relative power of 0.89 at 3 Hz during tACS steady-state entrainment period and hyperpolarizing-only stimulation, right, with average relative power of 0.98 at 3 Hz during tACS steady-state entrainment period).
Applying hyperpolarizing-only stimulation was more effective than depolarizing-only stimulation at entraining the network. A, Network activity for the PY layers of networks stimulated using depolarizing-only (left) and hyperpolarizing-only (right) stimulation waveforms. Both networks were entrained to the stimulation frequency, but the network stimulated with the hyperpolarizing-only waveform exhibited more regular activity. B, Two-dimensional maps of the PY layers of the two networks indicating membrane voltage (left), firing rate (middle), and synaptic depression (right) of each neuron. The network subjected to hyperpolarizing-only stimulation exhibited a high density of hotspots, which were able to coalesce more quickly than the depolarizing-only stimulated network.
We then compared the 2-D maps of membrane voltage, firing rate, and synaptic depression of both PY networks to further understand the role of spatiotemporal structure in the differential effect of depolarizing-only and hyperpolarizing-only tACS (Fig. 5B, 2-D maps of PY layer for networks receiving depolarizing-only, left, and hyperpolarizing-only, right, stimulation). In the network with depolarizing-only stimulation, two different activity modes occurred consecutively. First, just before the onset of the depolarizing phase of stimulation (simulation time = 1.3 s), one or two hotspots of activity formed in areas of low synaptic depression (high value of D). These hotspots slowly expanded outward, similar to the initial hotspots of the unstimulated or tDCS-stimulated networks. The second activity mode occurred after the onset of the depolarizing phase of stimulation, in which areas of the network not yet activated by the wavefront rapidly formed multiple hotspots that coalesced, similar to the behavior seen during full-wave tACS. In summary, the occurrence of “premature hotspots” before the onset of the depolarizing phase of stimulation disrupted the regularity of network oscillations.
For hyperpolarizing-only stimulation, we found exclusively this second, rapidly synchronizing, activity mode. Importantly, there were also areas of low synaptic depression at stimulation phase slightly before 2 π (i.e., just before the depolarizing phase of stimulation would occur with full-wave tACS). However, unlike the behavior of the depolarizing-only stimulated network, hotspots of activity were prevented from prematurely forming by the hyperpolarizing stimulation. As a result, a flood of uniformly distributed hotspots followed the offset of stimulation (similar to what we found stimulating with full-wave tACS). Therefore, this represents the underlying mechanism of hyperpolarization-induced entrainment by tACS.
Observing network and hotspot activity as a function of the phase of stimulation revealed further similarities and differences between the stimulation paradigms. While the network activity was uniformly distributed across the phase of stimulation for the unstimulated network, both full-wave tACS and hyperpolarizing-only stimulated networks had highly regular network activity occurring almost exclusively during stimulation phases ranging from 0 to π (i.e., when depolarizing, or no stimulation was applied). Although the depolarizing-only stimulated network exhibited much more regularity than the unstimulated conditions, the PY activity was more variable than that of the full-wave tACS and hyperpolarizing-only networks (Fig. 6A, network activity as a function of stimulation phase; left, unstimulated; left-middle, full-wave tACS; right-middle, depolarizing-only; right, hyperpolarizing-only). Additionally, PY activity for the depolarizing-only stimulated network continued during what was the DOWN state in the other two stimulated networks (in particular, for phases from 1.5 to 2 π, PYs were active 31% of the time, compared with only 3.2% and 5.6% of the time for tACS and hyperpolarizing stimulated networks). We also found that for both hyperpolarizing-only and full-wave tACS stimulated networks, the number of hotspots peaked at both onset and offset of stimulation, similar to the bimodal peaks observed in the time course of hotspot activity for tACS stimulated network (Fig. 6B, hotspot count as a function of stimulation phase). For the depolarizing-only stimulated network, the distribution of hotspots was quite uniform, with a slight peak during offset of network activity close to the stimulation phase of π. As expected from the network activity, we found that the majority of onset sites (defined as hotspots occurring while network activity was <10% of peak network activity) for the other two stimulated networks occurred close to a stimulation phase of 0 (entropy of both distributions = 0.07, whereas the distributions of the unstimulated and depolarized-only networks had entropies of 0.50 and 0.12, respectively; Fig. 6C, onset sites as a function of stimulation phase).
Phase space plots of networks, which were unstimulated (left), received full-wave tACS (middle left), administered depolarizing-only stimulation (middle right), and hyperpolarizing-only stimulation (right). Networks receiving full-wave tACS and hyperpolarizing-only stimulation are nearly identical, demonstrating that these stimulation paradigms had equivalent effects on the network. A, Network activity of the PY layer as a function of the stimulation phase. PYs were active almost exclusively during phases 0 to π of the stimulation for hyperpolarizing-only stimulation and full-wave tACS networks and was at a maximum when the phase of stimulation was ∼0.5 π. PY activity of the depolarizing-only stimulated network was more widely distributed and less regular. B, Hotspot count as a function of the stimulation phase. C, Histogram of total number of onset sites. Onset sites were defined as hotspots occurring when network activity was increasing but <10% of maximum network activity.
Further analysis of the network parameters revealed that the identical results observed for the hyperpolarization-only and full-wave tACS stimulated networks were due to the regulation of synaptic depression by the stimulation. We found that the relationship between synaptic depression and PY activity was tightly coupled for only these two stimulation paradigms (Fig. 7A, average synaptic depression across all PYs as a function of percentage PYs firing). For the hyperpolarization-only and full-wave tACS stimulated networks, the onset of stimulation and the consequential entrainment of the network resulted in an increase in the local maxima of the network's synaptic depression variable, D (Fig. 7B, local maxima of average synaptic depression variable of all PYs across time; median values during the period of steady-state entrainment with tACS were 0.36 ± 0.06, 0.478 ± 0.007, 0.43 ± 0.02, and 0.468 ± 0.008 for unstimulated, tACS, depolarization-only, and hyperpolarization-only stimulated networks, respectively).
Correlation of synaptic depression and network activity. A, The relationship between average synaptic depression of PYs to the network activity of PYs, for an unstimulated network, a network stimulated using full-wave tACS, depolarizing-only stimulated network, and hyperpolarizing-only stimulated network (top to bottom). Networks with more regular cortical dynamics (full-wave tACS and hyperpolarizing-only stimulation) experienced tighter coupling between the two parameters. B, Local maxima of average synaptic depression throughout simulation. Stimulation onset occurred at time 0. As the network became more regular, synaptic depression was able to reset more fully before it was depressed by the next UP state.
Sparse global stimulation was more effective than spatially localized stimulation at entraining the network
To validate the robustness of these findings, we simulated additional stimulation conditions that may occur when applying tACS. First, we included the stimulation of INs in addition to PYs and found that the results were nearly identical to stimulation of just PYs (relative power during stimulation 0.998 vs 1 for only PYs and both PY and IN stimulation paradigms, respectively). Next, we applied tACS and stimulated only the PYs as before; however, now the magnitude of stimulation received by each neuron was drawn from a uniform distribution with the same mean value as the stimulation amplitude. As the bounds of the distribution were increased, we found that the response of the network remained the same (relative power during stimulation remained at 1 as bounds of distribution increased from ±10% to ±100% of stimulation amplitude). Thus, neither targeting both cell types nor introducing heterogeneity in the stimulation amplitude across PYs significantly affect the results shown so far.
We next applied a spatially localized (SL) tACS paradigm in which only a percentage of the PYs were stimulated (Fig. 8A, inset, 2-D map of PYs for a 50% SL stimulated network). Here, we found that neurons receiving stimulation successfully entrained to the stimulation frequency; however, the unstimulated neurons were out of phase with the stimulation waveform (Fig. 8A; overall activity for SL stimulated network with 50% of PYs stimulated; stimulated neuron activity is shown in dark blue, and unstimulated neuron activity is shown in light blue). Closer inspection revealed that this differential effect of SL stimulation on the stimulated and unstimulated neurons was due to the spatiotemporal dynamics of the network (Fig. 8B, 2-D maps indicating the firing rate of individual PYs of the network from Fig. 8A at different times during an UP state; the stimulated top half and unstimulated bottom half of the network are separated by a red line). After steady-state entrainment was reached in the area of stimulated neurons, UP state initiation was characterized by a high density of hotspots in the area of stimulated neurons, which quickly converged and created regular high-amplitude UP states. However, in the unstimulated area, UP states resulted from a slow propagating wave. Therefore, the phase offset of the unstimulated neurons relative to the stimulation waveform was caused by the tight local coupling of the network. Specifically, entrainment of the unstimulated part of the network was mediated by propagating waves that originated in the stimulated part of the network.
Comparison of SL and RG stimulation. A, Overall PY activity for a SL stimulated network with 50% of PYs stimulated. Stimulated neuron activity (dark blue) and unstimulated neuron activity (light blue) are shown. Inset, A 2-D map of PYs with stimulated neurons shown in black. B, 2-D maps of PYs of SL stimulated network with 50% of PYs stimulated indicating firing rate of individual neurons from PYs of network from A at different times through an UP state. The stimulated (top) and unstimulated (bottom) halves of the network are separated by a red line. C, Total activity for RG stimulated network with 50% of PYs stimulated. D, 2-D maps of the PY subnetwork of the RG-stimulated network indicating the firing rate of each neuron with each frame taken at a different time during the initiation of an UP state. E, The relative power of the network activity versus the percentage of PYs stimulated for networks receiving SL (left) and RG (right) stimulation.
We next studied a random-global (RG) stimulation paradigm in which the stimulated PYs were distributed randomly throughout the layer. This stimulation modality resulted in identical activity for both stimulated and unstimulated neurons (Fig. 8C, total activity for RG stimulated network with 50% of PYs stimulated). Due to the global distribution of stimulated neurons and the tight local connectivity of the network, it was likely that an unstimulated neuron had neighbors that were stimulated. Therefore, the stimulation had a more direct effect on unstimulated neurons than in the SL paradigm. This resulted in the formation of a high density of hotspots throughout the entire PY layer during the onset of an UP state (Fig. 8D, 2-D maps indicating the firing rate of each PY in the RG stimulated network; each frame was obtained at a different time during the initiation of an UP state).
We quantified the effectiveness of SL versus RG stimulation by conducting an additional set of simulations across eight different networks. We applied both SL and RG stimulation to each of the eight networks to parameterize the percentage of neurons stimulated using each modality and quantified the power of total PY activity during the period of stimulation (Fig. 8E, the power of PY activity versus the percentage of PYs stimulated for networks receiving SL, left, and RG, right, stimulation). We found that for SL stimulation, the ability to entrain a network by stimulating only a subset of neurons varied across networks. It was only when all neurons were stimulated that the networks were reliably entrained (SD of power across networks at 10% = ±0.283, at 50% = ±0.181, and at 90% = ±0.160). However, we found that RG stimulation was more robust and all networks reliably entrained when ≥50% of the PYs were stimulated (SD of power across networks at 10% = ±0.080, at 50% = ±0.033, and at 90% = ±0.026). Therefore, the mechanisms of entrainment are robust to the details of the stimulation paradigm.
Entrainment was mediated by resonance dynamics
We next performed a comprehensive array of simulations to evaluate the degree of entrainment of the endogenous oscillation by tACS as a function of stimulation amplitude and frequency (Fig. 9A, tACS applied as a function of amplitude, from 1 to 13 pA, and frequency, from 0 to 6 Hz; the relative power at stimulation frequency is color coded). Dynamic systems theory predicts that systems with intrinsic periodic dynamics have preferred stimulation frequencies at which weak periodic perturbations are particularly effective at entraining the system (termed “resonance frequencies”). Indeed, tACS with matched frequency at 3 Hz was most effective in entraining the network at the lowest amplitude (average relative power of 0.97 at amplitude of 5 pA during period of entrainment); the harmonic frequency 6 Hz represented a second stimulation paradigm that also easily entrained the network (average relative power of 0.79 at amplitude of 5 pA during the period of entrainment). Other stimulation frequencies had less of an effect on network oscillations when applied at weak amplitudes (the mean of the average relative power of 0.47 at an amplitude of 5 pA during the period of entrainment for stimulation at all frequencies excluding 3 and 6 Hz). At higher stimulation amplitudes, the range of stimulation frequencies that successfully entrained the network broadened (Fig. 9B, color-coded frequency that exhibited maximum power during tACS). This relationship between the frequency of a periodic perturbation and its effect on an intrinsic oscillator is an Arnold tongue, a frequently observed phenomenon in dynamic systems (Pikovsky et al., 2001). We also found that stimulation paradigms that successfully entrained the network had a higher average number of hotspots (Fig. 9C, average number of hotspots is color coded), indicating that the above-described mechanism of entrainment was conserved for all effective tACS waveforms.
Series of simulations with tACS applied as a function of stimulation amplitude (1–13 pA) and frequency (0–6 Hz). A, Network activity of the PY layer for each simulation. Warmer colors indicate the higher power of network activity at the stimulated frequency. A prominent Arnold's tongue was present, in which entrainment of the network at its intrinsic resonant frequencies (3 and 6 Hz) required the lowest amplitude of stimulation. The range of frequencies the network was able to entrain increased with increasing stimulation amplitude. B, Frequency of maximum power for the stimulated networks. Same structure as in A. C, Average number of hotspots for the stimulated networks.
To quantify the entrainment stability, we conducted longer simulations with 25 s of stimulation (Fig. 10, spectrograms, with stimulation frequencies ranging from 2 to 6 Hz and relative power color coded). For the stimulation amplitude and frequency pairs that were within the Arnold tongue, entrainment was stable over time (SD of 0.014 for power at 3 Hz during the period of entrainment with 3 Hz tACS at 9 pA amplitude). In all other cases, we found slow patterning between epochs of stable entrainment and desynchronized activity. The relative time spent entrained versus nonentrained increased with (1) stimulation amplitude and (2) decreasing difference between stimulation and intrinsic oscillation frequency.
The presence of network bistability with alternating periods of entrainment and lack of entrainment for stimulation frequencies that do not match intrinsic (harmonic) frequencies. The same matrix representation as in Figure 9, with spectrograms shown for each stimulation condition. Fragmentation of power away from intrinsic frequencies resulted in macroscopic, bistable dynamics with periods of entrainment interleaved with periods of seemingly little stimulation effect.
tACS enhanced cortical oscillation in vivo in a frequency-specific manner
So far, we have presented findings from large-scale computer simulations that demonstrated the capability of tACS to enhance endogenous cortical oscillations with frequencies close to the stimulation frequency. We next asked whether tACS has the same effect in vivo. To this purpose, we performed multichannel electrophysiological measurements of action potential firing (multiunit activity) in anesthetized ferrets (N = 2 animals, 480 trials per stimulation frequency). In contrast to rodents, ferrets exhibit a convoluted cortex, which makes the effects of TCS less predictable, as the orientation of the somato–dendritic axis of neurons relative to the weak electric field generated by TCS varies due to the topology of gyri and sulci. We anesthetized the animals to generate a slow, rhythmic baseline state in the absence of stimulation. tACS was applied through bilaterally positioned large-stimulation electrodes (Fig. 11A, top). The application of stimulation increased the rhythmic structure of the multiunit firing (Fig. 11A, bottom). We studied the frequency spectrum of the smoothed multiunit activity, allowing us to probe the instantaneous effect of stimulation, in contrast to human EEG studies in which changes in oscillation power can only be assessed after the end of stimulation due to electrical artifact issues. Averaged across trials (30 s control, 30 s stimulation, 30 s control), we found a pronounced increase in rhythmic power at the stimulation frequency (Fig. 11B, stimulation frequencies 0.5–2.5 Hz, left to right, averaged across recording channels on depth probe). When analyzed individually per recording site, we found that this spectral enhancement was mostly limited to superficial sites, corresponding to layers I–II/III (Fig. 10C, spectrograms as a function of depth, stimulation frequencies from 0.5 to 2.5 Hz, top to bottom). There was no pronounced outlasting effect after stimulation ended. For the 480 trials at each stimulation frequency, we calculated ratios of spontaneous activity after stimulation to before stimulation (0.5 Hz stimulation ratio = 1.0319, p = 0.010; 1.0 Hz stimulation ratio = 1.0057, p = 0.2972; 1.5 Hz stimulation ratio = 0.9976, p = 0.5793; 2.0 Hz stimulation ratio = 0.9915, p = 0.0062; 2.5 Hz stimulation ratio = 1.0038, p = 0.4612). The time-averaged spectra (Fig. 11D, fs = 0.5–2.5 Hz, top to bottom) further confirmed localized power enhancement at the stimulation frequency (median in black, SEM in gray). Interestingly, for slow fs values (0.5 and 1.0 Hz), the effect of stimulation appeared to be more broadly distributed, potentially indicating that these frequencies did not match the intrinsic frequency. Across all trials, all stimulation conditions exhibited significant enhancement of the power at the stimulation frequency (Fig. 11E, scatter plot of all trials, fs is color coded; median enhancement for fs = 0.5 to 2.5 Hz at superficial electrode sites 1.11, 1.15, 1.15, 1.14, and 1.13, and at deep electrode sites 1.08, 1.06, 1.04, 1.08, and 1.12, all p < 0.05, Wilcoxon rank-sum test). Not only did tACS enhance the oscillation apparent in the multiunit firing, but tACS also phase locked multiunit firing to the stimulation sine wave (Fig. 11F, phase histograms as a function of depth and fs, entropy is color coded and hotter colors indicate lower entropy and therefore more pronounced effect of tACS). Together, these experimental findings demonstrate that, indeed, very weak and global electric fields have similar overall effects to our computational model in complex gyrencephalic brains.
tACS in anesthetized ferrets enhanced cortical oscillations at the stimulation frequency.A, Top, Schematic representation of experimental setup; bottom, sample multiunit activity and stimulation trace. B, Averaged spectrogram for all stimulation frequencies. C, Depth profiles of spectrograms. D, Time-averaged change in spectra (gray: 95% confidence interval). E, Power of network activity at stimulation frequency during versus before tACS. Stimulation frequencies are color coded. F, Phase histograms show the preferred phase of network activity for superficial cortical layers (entropy is color coded, with hotter colors indicating lower entropy.).
Discussion
Recent tACS studies have shown promising effects for the targeted modulation of cortical oscillations and associated motor and cognitive functions (Antal et al., 2008; Kanai et al., 2008, 2010; Moliadze et al., 2010; Zaghi et al., 2010; Feurra et al., 2011a; Laczó et al., 2012). In this study, we used large-scale computational simulations to elucidate the mechanisms by which global weak perturbations such as TCS can modulate macroscopic network dynamics. The size of the computational model (total of 200,000 neurons) enabled us to study synchronization in a network that intrinsically oscillated at ∼3 Hz and exhibited propagating activity across the PY network. Our results highlight the benefits of using such large-scale models since the main mechanism of modulating the number of near-simultaneously occurring hotspots depends on a large neuronal territory where such events can occur. Interestingly, our results suggest that the overall interaction dynamics of the network with the applied tACS can be reduced to a very simple oscillator driven by a periodic force (Pikovsky et al., 2001). Evidence for this is provided in a recent landmark study by Polanía et al. (2012) where in-phase tACS at 6 Hz induced frontoparietal theta synchronization and enhanced working memory performance, while out-of-phase 6 Hz tACS decreased performance. This causal link between theta-band synchronization and working memory was further confirmed by the absence of effects of 35 Hz tACS. Thus, this study demonstrated phase- and frequency-specific effects of tACS in agreement with our findings on resonance dynamics.
Although our study exclusively focused on slow cortical rhythms, it can be argued that our findings on the underlying mechanisms of weak global perturbations are more broadly applicable. In essence, TCS enhanced the endogenous oscillation in our simulations by controlling the number of hotspots or initiation sites that then rapidly expanded by lateral propagation through local excitatory connections. Several locations in the network were very close to the threshold for massive positive feedback that produced hotspots due to the local PY–PY coupling. TCS then pushed these extra “potential hotspots” above the threshold for becoming full-blown initiation sites. Similar thresholds and nonlinearities commonly exist in oscillating neuronal networks in which epochs of quiescence alternate with epochs of neuronal firing (e.g., cortical “slow oscillation”). Therefore, the dynamics identified in this study may well generalize to different network configurations and endogenous oscillation frequencies. In this study, the model networks exhibited pronounced sensitivity to weak perturbations. This mediated the entrainment of the intrinsic oscillations by tACS. Specifically, our results suggest that increasing the activity of globally distributed and highly connected neurons is sufficient to trigger a transition between a DOWN and an UP state, two very different activity modes. Interestingly, several recent reports of electrophysiological studies in vivo point toward such high sensitivity of networks to perturbations of just a single cell (Li et al., 2009; London et al., 2010). These results bolster the biological plausibility of such high sensitivity to weak global perturbations.
We found that short-term synaptic depression enabled TCS to exert pronounced effects on cortical oscillations. In particular, the hyperpolarizing phase of the stimulation sine wave kept the network quiet such that synaptic depression of all excitatory synapses could recover, allowing for a higher-gain positive-feedback PY–PY excitatory interaction for the next UP state (and therefore enhanced synchronization). The presence and strength of short-term synaptic depression in cortex in vivo is debated and may be different across physiological conditions (Zucker and Regehr, 2002; Reig et al., 2006). However, any other negative-feedback mechanism that has a comparable recovery time constant could play the same role in the dynamics. Therefore, our findings are not exclusively valid for scenarios in which excitatory synapse strength depresses because of repeated activation. For example, BK potassium channels (Lee and Cui, 2010) that are gated by intracellular calcium (Ca2+) concentration, which increases during UP states, can display similar dynamics role since keeping activity low via stimulation would allow increased Ca2+ extrusion and result in the closing of BK channels (Fröhlich and Bazhenov, 2006). In summary, several other activity-dependent modulators of intrinsic or synaptic excitability are expected to play a very similar role in mediating the effects of tACS on endogenous network oscillations.
Our results demonstrate the higher efficacy of tACS versus tDCS for the instantaneous entrainment of endogenous cortical oscillations. Our model does not include long-term potentiation or other plasticity mechanisms with longer time-scales and can therefore not address the mechanisms of the ongoing effects of TCS following the offset of stimulation. Furthermore, our results on resonance dynamics provide a note of caution for the choice of stimulation frequency since cortical oscillations were strongly suppressed for stimulation frequencies in between the dominant endogenous frequency and its first harmonic frequency. An important implication of this finding is that the frequency of applied stimulation should be matched to the frequency of the endogenous oscillatory state. However, macroscopic activity such as that measured by EEG typically exhibits a 1/f spectrum without any single strong peak (Allegrini et al., 2009) such as in our model. Therefore, the choice of stimulation frequency could represent a serious challenge as there is no clear preferred resonance or peak frequency. Beyond these technical concerns, our findings may advocate for closed-loop feedback systems similar to those already implemented in animal studies, in which the stimulation waveform is dynamically adjusted to ongoing activity patterns (Berényi et al., 2012).
Footnotes
This research was supported by the University of North Carolina (UNC) Department of Psychiatry, the UNC School of Medicine, and the Foundation of Hope. The authors thank Stephen Groner for help with the initial code development; Davis Bennett for help with the manuscript; Dr. Jan Prins for expert advice on parallel computing; and the UNC Scientific Computing group, in particular Dr. Mark Reed, for support with GPU-based high-performance computing.
The authors declare no competing financial interests.
- Correspondence should be addressed to Flavio Frohlich, 115 Mason Farm Road, NRB 4109F, Chapel Hill, NC 27599. flavio_frohlich{at}med.unc.edu
