Abstract
Rhythmic brain activity plays an important role in neural processing and behavior. Features of these oscillations, including amplitude, phase, and spectrum, can be influenced by internal states (e.g., shifts in arousal, attention or cognitive ability) or external stimulation. Electromagnetic stimulation techniques such as transcranial magnetic stimulation, transcranial direct current stimulation, and transcranial alternating current stimulation are used increasingly in both research and clinical settings. Currently, the mechanisms whereby time-dependent external stimuli influence population-scale oscillations remain poorly understood. Here, we provide computational insights regarding the mapping between periodic pulsatile stimulation parameters such as amplitude and frequency and the response dynamics of recurrent, nonlinear spiking neural networks. Using a cortical model built of excitatory and inhibitory neurons, we explored a wide range of stimulation intensities and frequencies systematically. Our results suggest that rhythmic stimulation can form the basis of a control paradigm in which one can manipulate the intrinsic oscillatory properties of driven networks via a plurality of input-driven mechanisms. Our results show that, in addition to resonance and entrainment, nonlinear acceleration is involved in shaping the rhythmic response of our modeled network. Such nonlinear acceleration of spontaneous and synchronous oscillatory activity in a neural network occurs in regimes of intense, high-frequency rhythmic stimulation. These results open new perspectives on the manipulation of synchronous neural activity for basic and clinical research.
SIGNIFICANCE STATEMENT Oscillatory activity is widely recognized as a core mechanism for information transmission within and between brain circuits. Noninvasive stimulation methods can shape this activity, something that is increasingly capitalized upon in basic research and clinical practice. Here, we provide computational insights on the mechanistic bases for such effects. Our results show that rhythmic stimulation forms the basis of a control paradigm in which one can manipulate the intrinsic oscillatory properties of driven networks via a plurality of input-driven mechanisms. In addition to resonance and entrainment, nonlinear acceleration is involved in shaping the rhythmic response of our modeled network, particularly in regimes of high-frequency rhythmic stimulation. These results open new perspectives on the manipulation of synchronous neural activity for basic and clinical research.
Introduction
Rhythmic neural dynamics are a mechanism through which internal states exercise (top-down) influences on stimulus processing and in turn impact perception (Engel et al., 2001; Varela et al., 2001). Synchronization is considered one such mechanism, operating across multiple scales of brain circuitry from single units to interregional coherence (Singer and Gray, 1995; Engel and Singer, 2001). The functional role of oscillations within the gamma band (30 Hz–100 Hz) has received considerable attention (Whittington et al., 1995; Wang and Buzsáki, 1996; Jadi and Sejnowski, 2014). More recently, lower frequencies, which likely rely on different cellular mechanisms and operate over longer timescales than their high-frequency counterparts (Kajikawa and Schroeder, 2015; Podvalny et al., 2015), have been the focus of increasingly intensive research (Lorinscz et al., 2009; Haegens et al., 2014). In particular, the alpha frequency (8–12 Hz) is known to correlate strongly with cognitive performance and serves as a robust predictor of the degree of cognitive decline between healthy and pathological groups (Klimesch, 1999; Başar, 2012; Basar and Guntekin, 2012; Vlahou et al., 2014). The alpha peak frequency has also been found to predict the temporal resolution of visual perception (Samaha and Postle, 2015). Therefore, numerous studies have investigated the consequences of variations in alpha oscillations on cognition, attention, and perception to increase our understanding of the reciprocal interactions between slow oscillatory activity and population neural coding (Lakatos et al., 2008; Siegel et al., 2012). For example, shifts in the peak alpha frequency of different subjects (also called individual alpha frequency) are observed reliably when attentional demands increase during visual tasks (Samaha et al., 2015), overt behavior (Gutmann et al., 2015), and/or sensorimotor inputs (Hülsdünker et al., 2015).
These observations have raised the prospect of manipulating (alpha) rhythmic activity in a controlled manner either through attention-related processes or brain stimulation methods (Thut et al., 2012; Frohlich, 2015). It has indeed been demonstrated that transcranial magnetic stimulation (TMS) and transcranial direct current stimulation and transcranial alternating current stimulation (tACS) can perturb individual alpha oscillations and have a direct effect on visual stimulus perception (Chanes et al., 2013; Cecere et al., 2015) and task performance by reinforcing endogenous slow-wave rhythms (Klimesch, 1999; Henry et al., 2014).
One implication of such findings is that brain stimulation can result in the selective engagement of endogenous (intrinsic) oscillations, effectively making it a potential means of manipulating and controlling cognition and treating neurobiological disorders. This critical question has received much attention in recent studies (Frohlich, 2015; Witkowski et al., 2015) in which mechanisms such as resonance (Thut et al., 2012; Ali et al., 2013) and entrainment (Thut et al., 2011) have been implicated and found to depend strongly on the ongoing dynamics of the stimulated networks (Neuling et al., 2013; Alagapan et al., 2016).
In order for brain stimulation approaches to be optimally effective, it is necessary to have a better mechanistic understanding of their consequences on brain activity, something that can be garnered from computational modeling. However, computational modeling of this issue has rarely considered a wide range of stimulation parameters (Reato et al., 2010). Due to the nonlinearity of neural circuits and the large number of possible stimulation patterns, the mechanisms involved in the combination of exogenous and endogenous oscillatory signals are furthermore poorly understood. To address these limitations, we here consider the response of a simplified cortical network to periodic pulsatile stimuli. Specifically, using a generic computational model of spiking excitatory and inhibitory neurons that exhibits spontaneous alpha oscillatory activity (here set at 10 Hz), we explore different combinations of stimulation frequencies and amplitudes systematically, revealing those in which resonance and/or rhythmic entrainment can be evoked. Then, building on recent advances in neurodynamics (Lefebvre et al., 2015), we report novel regimes of oscillatory acceleration that characterize network responses to input stimuli of high frequencies (50–100 Hz). Our results provide new perspectives on the response of synchronous nonlinear neural systems, in which a plurality of linear and nonlinear mechanisms are combined to shape intrinsic alpha-like oscillatory activity.
Materials and Methods
Cortical network.
To highlight the influence of stimulation on local oscillatory dynamics, we deliberately chose a relatively well understood hybrid Amari model (Amari, 1977) of a cortical network and evaluated its behavior when subjected to dynamic stimulation, that is, during forcing. Our model takes into account the variability of neuronal responses, yet is simplistic enough to remain both computationally and analytically tractable. It consists of a network of pyramidal excitatory neurons (e) and inhibitory interneurons (i), the individual spiking activity of which follows the nonhomogeneous Poisson processes as follows:
where Xe,ij(t) = ∑{tl}δe,ij(t − tl) is the spike train of the jth neuron in the population e or i. For simplicity, the activation function is defined by f[u] = 1 whenever u > h and 0 otherwise. It represents how the firing rate of the neurons relates to the cellular potential, u. The excitatory and inhibitory somatic membrane potentials uej(t) and uij(t) evolve according to the set of scalar nonlinear stochastic equations as follows:
Here, Ne and Ni refer to the number of e versus i cells, which follow the 4:1 ratio generally observed in cortical tissue (Koch, 1999). The rates αeu, αiu define the timescale of the respective somatic membrane. Adaptation currents vej(t) and vij(t) with gain b and rate a are also included to represent the habituation of neurons to both steady and fluctuating inputs (Prescott and Sejnowski, 2008).
Both excitatory and inhibitory populations are driven by an external stimulation signal I(t), which consists of pulsatile stimuli. In this model, recurrent inputs to the neurons take the form of postsynaptic potentials due to incoming spikes. The cross-population recurrent inputs Gnmj(t) are defined by the following:
where EPSPk(t) and IPSPk(t) are afferent postsynaptic excitatory and inhibitory potentials, respectively. They are computed by convoluting spike trains with exponential synapses with time constant τm as follows:
Both populations are further subjected to endogenous sources of noise ϕe,ij(t), which is assumed to follow spatially and temporally independent Gaussian white noise profiles with fixed amplitude D. Synaptic weights within (Weejk(c), Wiijk(c)) and between (Weijk(c), Wiejk(c)) excitatory and inhibitory populations exhibit sparse exponential profiles (Hellwig, 2000) with connection probability c as follows:
where σn,m2 = σe2, σi2 is the range of the excitatory and/or inhibitory interactions and x(k) refers to the spatial location of neuron k in the network. To represent the sparseness of cortical connections, synaptic weights wnmo(c) are set to 0 with probability of 1 − c. The matrix of synaptic weights is plotted in Figure 1A. The network extends in the one-dimensional spatial domain Ω and the neural lattice is built by distributing cells randomly according to a uniform distribution. Given the spatial organization of the network, synaptic interactions between neurons k and l are thus subjected to propagation delays τjk = |x(k) − x(j)|v−1, with v being the axonal conduction velocity, set here to 0.15 m/s. This value for the conduction velocity is realistic for intracortical connections (Hutt et al., 2003). This gives rise to a distribution of propagation latencies plotted in Figure 1B. For the parameters selected, neurons engage spontaneously in synchronized oscillatory alpha activity at a frequency of ∼10 Hz. Sample spiking responses for randomly selected inhibitory and excitatory cells are plotted in Figure 1C, alongside membrane potential traces exhibiting subthreshold oscillations within the alpha band. These oscillations result from the combined influence of propagation delays and slow cellular adaptation, which is known to lead recurrent nonlinear networks in a state of robust synchrony (Lefebvre et al., 2011). This can be seen from the full network activity plotted in Figure 1D, where the spiking of the neurons is modulated by a global rhythm of 10 Hz. Model parameters are summarized in Table 1 below. We note that, throughout, the network parameters were left unchanged and only stimulation settings such as input frequency and amplitude were varied. Being a single-compartment, conductance-based model, the system above shares many similarities with the well known two-population integrate-and-fire networks that have also been shown to sustain global synchronous activity (Brunel, 2000; Buzsáki and Wang, 2012). Key differences lie in the spike initiation mechanism, which in our case adds an additional source of variability in the firing of the modeled neurons, and in the presence of slow inhibition, resulting in slower oscillations frequencies.
Network connectivity and intrinsic activity without stimulation. A, Matrix of synaptic weights that define the connectivity of the network within (Wee, Wii) and between (Wie, Wei) excitatory and inhibitory neural populations. As defined, the connectivity is spatially profiled, where proximal (resp. distal) neurons share stronger (resp. weaker) synaptic connections. B, Distribution of propagation delays in the network. Given a finite conduction velocity of 0.15 m/s, the timing of synaptic interactions is delayed. The distribution shows that the network dynamics are dominated by many short latencies and few longer ones. C, Sample activity of randomly selected excitatory (black) and inhibitory (blue) neurons. The spiking activities Xej and Xik of these neurons are shown above the corresponding membrane fluctuations, which exhibit clear, albeit noisy, subthreshold oscillations. D, Network intrinsic spiking activity without stimulation. Both populations are set in stable synchronous activity with firing rates modulations of ∼10 Hz.
Model parameters
Spectral analysis.
Spectral analysis was performed using a fast Fourier transform routine using freely available C++ scripts (Press et al., 2007).
Stimulation.
The stimulation I(t) in Equation 1 is a train of phasic pulses of duration δt = 300 μs of given forcing frequency ωf and intensity If. In the analysis, both frequency and intensity are changed to explore a variety of independent combinations. This stimulus offers a richer spectral signature compared with sinusoidal signals, for which only one frequency is represented. Later on, it is useful to compute the time average of the stimulation pattern μI = <I(t) > T over a given epoch that is given by the following:
which also allows also the computation of the variance as Var(I) = <(I(t) − μI)2 > T = IfμI − μI2. We deliberately chose short stimulation epochs (i.e., 1 s) compared with the timescale usually involved during continuous stimulation paradigms in clinical settings (lasting several minutes) because we were investigating the online effect during stimulation, which has been shown to exist for brief stimulation duration of a few seconds (Fröhlich and McCormick, 2010). For the investigation of plastic effects, such as long-term synaptic potentiation (Huerta and Volpe, 2009), longer stimulation duration would be required.
Neuroelectric activity.
In our model, the network neuroelectric activity is determined by a weighted sum over somatic excitatory and inhibitory potentials as follows:
where φe,ik are real positive coefficients. Here, we assume that the network fine scale structure is unknown, and thus consider random weights as follows: φe,ik = [0,1]. This is done to take into account, in our calculations, of different sources of observational variability and constraints that we do not model explicitly while avoiding any filtering issues in the resulting spectral calculations. However, specific choices of the φe,ik distributions can be made to increase the similarity of the network activity to physiological signals such as LFPs and EEG (Lindén et al., 2010; Mazzoni et al., 2015).
Spike coherence.
To compute spike coherence between spike trains (Gabbianni and Koch, 1998), we selected randomly two excitatory neurons, and j, within the cortical populations and binned their respective responses over a given time window ΔT. Therefore, for neuron j, Xjbin(kΔT) = 1 if a spike occurred during the interval [kΔT, (k + 1)ΔT] and zero otherwise, where k = T/ΔT. Coherence, θ(ΔT), was then computed as follows:
where ∑T denotes the sum over time intervals. For θ(ΔT) = 0, two spike trains are independent and θ(ΔT) = 1 reflects maximum coherence. A nonlinear (resp. linear) dependence of the coherence on the time window will appear whenever synchronous (resp. asynchronous) spiking exists in the data (Wang and Buzsaki, 1996; Hutt and Buhry, 2014). For the purposes of the current analysis, we have chosen ΔT = 10 ms. This value was selected to optimize the coherence difference between original (synchronous) and shuffled spike trains.
Mean coupling strength.
To measure the relative level of recurrence and nonlinearity in the system, we computed the network-average recurrent inputs received by neurons during stimulation epochs and how they changed with respect to different stimulation parameters. The mean coupling strength, referred here as the κ index, is defined by the following:
where ≪.≫N,T stands for the mean taken over the population of size N and during a time interval of duration T. For the purposes of the analysis, we have chosen to compute the mean coupling strength over the full duration of the stimulation epochs to average out time fluctuations conveyed by the oscillatory nature of the network output.
Results
To understand the response evoked during stimulation, we investigated the behavior of a generic cortical network expressing alpha synchrony and subjected to periodic stimulation of various combinations of frequencies and amplitudes and observed the evolution of the spectral properties exposed by the network's neuroelectric output. The goal here was not to derive a precise mathematical formulation of the resulting dynamics, but instead to provide a general overview of driven oscillatory responses unfolded in stimulation parameter space. Limited by the dimensionality of possible stimulation waveforms, we focused our attention on repetitive pulsatile stimuli that closely resemble transcranial magnetic stimulation waveforms and/or direct cortical microstimulation paradigms (see Materials and Methods). The stimulation was assumed to be global: all neurons in the network were subjected to the same driving pattern. The aim of the study was to extract the mapping that links input to output spectral signatures. Therefore, we assumed stimulation epochs to be short compared with synaptic potentiation and/or depression such that plasticity had no significant impact on the dynamics. Therefore, our findings characterize instantaneous responses arising due to the recurrent architecture of the network and the rhythmic nature of its ongoing state.
Figure 2 details different response mechanisms that can be evoked during pulsatile forcing in our model network. Without stimulation, the network is found in a synchronized state, exhibiting strongly correlated activity where firing rate fluctuations oscillate at an intrinsic—also called “natural”—frequency of ωo ≅ 10 Hz. Different stimulation forcing amplitudes (If) and frequencies (ωf) lead the network through different response regimes such as resonance, in which intrinsic oscillations are enhanced (Fig. 2B); entrainment, where the systems' dynamics are locked to the stimulation (Fig. 2C); and nonlinear acceleration (NLA), in which the stimulation causes a shift in the system's intrinsic frequency (Fig. 2D).
Responses of synchronized neural population to periodic pulsatile stimulation. The network neuroelectric output A(t) is shown on the left. The associated power spectral density (PSD) is plotted in the middle next to the excitatory neurons spiking activity, which are shown on the right. A, According to the set of selected parameters, the network engages spontaneously into synchronous alpha activity. In absence of stimulation, the endogenous rhythm stabilizes around an unperturbed intrinsic frequency of ωo ≅ 10 Hz, which modulates the noisy spiking activity of the neurons. Once global stimulation is applied to the network, different mechanisms are recruited and explain the observed responses. B, If the stimulation frequency ωf is close to the intrinsic frequency, its harmonics, and/or subharmonics, resonance occurs. Although the network response frequency ω′ remains stable—that is, ω′ ≈ 
ωo—the associated spectral power increases significantly. This is also reflected in the spiking activity of driven neurons: firing rates increase into tightly synchronized bursts of spike discharges. C, Whenever the stimulation amplitude increases beyond a certain threshold, intrinsic ongoing activity in entrained by the stimulation drive and the network oscillates at a frequency equal to the stimulations, its harmonics, and/or subharmonics; that is, ω′ ≈ 
ωf. In this example, a strong rhythmic stimulation of ωf ≅ 17.5 Hz was applied and harmonic entrainment can be observed in the 1:2 Arnold tongue as the network peak frequency shifts to ω′ ≅ 8.8 Hz. D, For high-frequency stimulation, nonlinear interactions provide the means of accelerating ongoing synchronous activity. A spectral shift Δ(ωfIf), the magnitude of which depends on stimulus parameters, arises due to a change of the network intrinsic oscillatory properties. In this example, a stimulation of ωf ≅ 97 Hz can be seen to accelerate the baseline activity from ωo ≅ 10 Hz up to ω′ ≅ 15 Hz.
To disambiguate network responses and see how they are shaped by the system's nonlinearity, we analyzed the network response during fully independent trial epochs of 1 s duration. During each epoch, we computed the power spectral density of the network neuroelectric activity (See Materials and Methods). Varying the stimulation frequency alone (i.e., for a fixed stimulation intensity), we found that the different mechanisms detailed in Figure 2 were evoked sequentially, yet in a nontrivial manner. Figure 3A displays the network‘s neuroelectric power spectral density (see Materials and Methods) evolving as a function of increasing stimulation frequency while keeping the stimulation intensity constant. As detailed in Figure 3B, alongside the fundamental frequency peak seen on the diagonal, one can see both harmonic and subharmonic peaks. For lower frequencies, the stimulation was found to tune the peak frequency expressed by synchronous neurons via a sequence of subharmonic entrainments. For high-frequency stimulation, the network was found in a state of fast oscillations in which ongoing cyclic activity was not entrained to the stimulation, but instead was accelerated with respect to the baseline intrinsic frequency.
Power spectrum of the network neuroelectric output A(t) during stimulation of fixed intensity but variable frequency. A, Increasing the frequency of the stimulation (ωf) shifts the peak frequency (ω′) and alters the spectral profile of the network's activity. The stimulation fundamental frequency is seen aligned on the diagonal. In addition, both harmonics (Hm:nω) and subharmonic (SHn:mω) of the stimulation frequency are visible (above and below the diagonal, respectively). For low-frequency forcing, the peak frequency is close to the system's intrinsic frequency ωo. However, entrainment occurs as ωf increases further, where the network peak frequency successively tracks the stimulation frequency ωf and its subharmonics. Once ωf exceeds 50 Hz, the peak frequency stabilizes to an accelerated frequency. B, Schematic representation of the spectral patterns seen in A. Shifts in the network peak frequency (ω′) are plotted in red. Subharmonic resonance can be observed as indicated by the black lines below the diagonal. Here, If = 1.5. In Figure 4, one can see how these patterns unfold as stimulation intensity is also changed.
How do the peak frequency and power change when stimulation intensity is varied? To better understand how the spectral patterns seen in Figure 3 unfold, we systematically varied both If and ωf and identified in each case the network peak frequency (ω′) and peak power at that same frequency. In Figure 4, multiple spectral patterns can be seen in both the power and frequency plots, indicating a plurality of network responses. For lower stimulation frequencies, the system's responses were expectedly found to be shaped by entrainment and resonance. Indeed, as seen in Figure 4A where the power of the peak frequency is plotted, one can see characteristic resonance peaks (R) aligned at multiples of the network intrinsic frequency (ωo) where the stimulation triggered the amplification of the endogenous alpha rhythm. This effect spanned all stimulus intensities. Harmonic resonances can also be seen where peaks appear at integer multiples of ωo (e.g., nωo) with gradually decreasing amplitudes. Nonlinear effects can also be seen readily. Resonance peaks, which would be vertically aligned at every multiple of the systems' intrinsic frequency in the case of a linear oscillator, are here arched. The effect, which is more salient for higher harmonics, is due to a shift in network intrinsic frequency due to the action of stimulation. In addition, one subharmonic resonance peak can be observed at ∼5 Hz. The different resonances expressed by the network are labeled in Figure 4B.
Spectral mapping between input stimulation parameters and the network neuroelectric output A(t). Stimulation frequency (ωf) and amplitude (If) were varied systematically while the network response was analyzed. A, Spectral peak power expressed by the network as a function of stimulation parameters. Characteristic resonance peaks, vertically aligned, are present at the system's intrinsic frequency (ωo ≅ 10 Hz) and at every integer multiple of it. Subharmonic resonance can also be seen at 5 Hz. Through these resonances, the system's output power is magnified due to the proximity of the stimulation frequency to the intrinsic frequency. The peaks, however, gradually vanish as the stimulation frequency increases >50 Hz, leading to responses of relative weak power for all remaining stimulation conditions. B, Labeled resonance (R) curves seen in A, shown along the subharmonic resonance at 5 Hz (SR). C, Network peak frequency (ω′) as a function of stimulation parameters. Arnold tongues representing entrainment are aligned vertically and indicate regimes where the system's frequency is locked to the stimulation frequency. Harmonic entrainment can also be seen. Entrainment occurs mainly for frequencies close the system's intrinsic frequency and its harmonics and vanishes as the stimulation frequency increases. For higher frequencies and amplitudes, however, nonlinear acceleration of intrinsic oscillatory activity can be observed. Such acceleration occurs for higher stimulation frequencies, where ω′ is seen to transit gradually to frequencies in the beta band by the action of forcing alone. D, Labeled entrainment (E) and NLA regimes as seen in C. Note that Arnold tongues overlap each other. The horizontal dashed lines at If = 1.5 seen in all panels represent the specific example detailed in Figure 3 in which the frequency, but not the intensity, was varied (within the interval 0–50 Hz). One can then see how peak power and frequency are shaped as stimulation intensity is also varied.
Other patterns emerge in Figure 4C, where frequency of the peak frequency is plotted with respect to varying stimulation parameters. For weaker stimulation frequencies, characteristic entrainment (E) patterns, also called Arnold tongues (Jensen et al., 1983), dominate the system's response: the network oscillatory behavior is locked to the stimulation frequency; that is: ω′ ≅ ωf. Subharmonic entrainment also occurs as the stimulation frequency increases: the network is entrained to oscillate at frequencies that are subharmonics of the driving frequency. Although similar to resonance peaks, the mechanism involved in generating these Arnold tongues is fundamentally different. We also note that. whereas the entrainment effect scales with stimulation intensity, it is nonetheless limited to driving frequencies that are close to ωo and/or its harmonics (Fig. 4A,B) and vanishes at higher stimulation frequencies.
For stimulation frequencies >50 Hz, another mechanism emerged and was found to dominate the response of the network, the dynamics of which enter a state of increasing recurrence. In this regime, high-frequency pulses increase the network response frequency via an NLA of the endogenous rhythmic activity. The intrinsic network frequency is shifted, leading to a gradual transition from the alpha (8–12 Hz) to the beta (12–30 Hz) band. In contrast to high-intensity/low-frequency stimulation, where resonance and entrainment dominate, NLA occurs in response to stimuli that have more subtle effects on the individual firing rates of driven neurons. Despite high forcing frequencies and amplitudes, NLA triggers rhythmic responses expressing peak power levels comparable to those seen in the unperturbed state (Fig. 2A) and dominates the peak frequency mapping at high frequencies. Visible entrainment and NLA regimes are labeled in Figure 4D.
Close scrutiny of the network firing rates during stimulation epochs also revealed the influence of stimulation on the spiking activity of individual neurons. We computed the normalized mean firing rates of excitatory neurons by first computing the network mean spiking activity in every stimulation condition and then comparing it with resting-state values. As seen in Figure 5A, for weaker stimulation frequencies (i.e., <50 Hz), neural discharge rates are strongly driven by network resonances: neurons' spiking activities track stimulation-induced changes in network peak power. Therefore, firing rates are strongly tied to the stimulation signal. However, for higher stimulation frequencies (i.e., >50 Hz), notably where NLA occurs, the network's mean firing rates still increase despite the weak neuroelectric power expressed (cf. Fig. 4A).
Firing rates and correlations as a function of varying stimulation parameters. A, The rhythmic stimulation was found to modulate individual firing rates non-monotonically. Changes in stimulation parameters triggered increases in spiking activity that mirrored the resonance peaks pattern found in Figure 3A. For higher frequencies and/or amplitude, a regime of high spiking activity can be found despite the absence of significant spectral power found in the same region. The mean firing rate of an excitatory neuron, averaged over a stimulation epoch of 1 s, is plotted. B, Spike coherence was also found to reflect resonance-induced modulations. Resonance was thus recruited to generate responses with highly correlated firing. Spike coherence was found to be low for high-frequency and/or high-amplitude stimuli despite strong spiking activity.
Therefore, to disambiguate the contributions of firing rates, power, and peak frequency to the emergence of synchronous dynamics in presence of stimulation, we further computed the response spike coherence between neurons as stimulation parameters were varied systematically. As seen in Figure 5B, the shadow of the system's resonance peaks were also found, indicating that rhythmic stimulation triggers, not only high firing rates via resonance, but also strongly correlated states, which is consistent with reported links between expressed neural correlations and the magnitude of spiking activity (De La Rocha et al., 2007). In contrast, weak to no spike coherence was found between the neuron's spike trains for stimuli of high frequency (Fig. 5B).
This stimulus-driven phenomenon echoes previous results that revealed the nonlinear mapping between stochastic input statistics and the spectral features of nonlinear recurrent networks (Lefebvre et al., 2015). Building on the timescale separation that subsists amid slow endogenous activity and high-frequency fluctuations, synergistic interactions between stimulation-induced fluctuations and network nonlinearity support a gain control mechanism in which the network susceptibility is enhanced with respect to baseline. In this framework, the stimulation recruits the network's recurrent nonlinearity and distorts its oscillatory properties, an effect that further scales with stimulation intensity; that is, ω′ ≈ ωo + Δ(If, ωf). To support this perspective, we calculated the stimulation time-average, [μ]I (see Materials and Methods) as stimulation settings were varied. As can be seen from Figure 6A, the stimulation mean intensity expectedly increases with both stimulation amplitude and frequency: maximal recruitment of network recurrent nonlinearity should thus mirror the μI distribution. To verify this, we computed in Figure 6B the mean coupling strength index κ (see Materials and Methods) as stimulation parameters were varied systematically. The κ index gauges the average intensity of recurrent interactions across the network and scales with network-induced intrinsic activity, as opposed to stimulation-driven activity; therefore, it is an indirect measure of the system's effective recurrence. The mean coupling strength was found to be shaped by stimulation distinctively compared with the neuron's firing rates and/or spike coherence: the κ index distribution correlated strongly with the NLA pattern shown in Figure 4C and did not build on from resonance and entrainment patterns despite being modulated by those as well. More importantly, maximal recruitment of network recurrence was achieved alongside maximal stimulation intensity, indicating that nonlinear effects support the oscillatory shifts observed in this region of parameter space. Together, these results suggest that, in addition to resonance and entrainment, rhythmic stimulation can also be used to recruit the intrinsic nonlinearity of cortical networks to shape the oscillatory features displayed by synchronous populations.
A shift in recurrence mediates the nonlinear acceleration observed for high-frequency stimulation. A, Stimulation time-average μI (see Materials and Methods) as stimulation settings are varied. As expected, μI scales with stimulation frequency and amplitude. B, The mean coupling strength κ (see Materials and Methods) follows the same trend and scales with stimulation parameters. Maximal recruitment of recurrent interactions is achieved with peak stimulation mean intensity and thus underlies the nonlinear acceleration of intrinsic oscillations.
Discussion
There is growing evidence demonstrating that noninvasive neurostimulation engages neural networks by perturbing endogenous oscillatory brain activity. Purported mechanisms include the generation of bursts of time-locked gamma oscillations (Huerta and Volpe, 2009) or the disruption and/or entrainment of endogenous alpha cyclic dynamics (Thut et al., 2011) that lead to quantifiable alterations in cognitive performance (Romei et al., 2011) and perception (Chanes et al., 2013; Cecere et al., 2015). Moreover, the influence of stimulation has further been reported to be state dependent (Neuling et al., 2013). The current quantitative assessment of neurostimulation's impact is oftentimes limited to frequency ranges that are close to the system's natural frequency in which resonance plays a dominant role. Here, we have explored the response of a model network of spiking neurons to various combinations of forcing intensity and frequencies to characterize better the plurality of network responses evoked.
Our results confirm that the engagement of nonlinear spiking networks to low-frequency forcing (i.e., <50 Hz) is shaped by resonance and entrainment. In this regime, an enhancement of the amplitude of the intrinsic oscillation (resonance) can be observed reliably if the stimulation frequency is close to the intrinsic frequency or its harmonics and subharmonics. For increased stimulation amplitudes, even if the network is stimulated at frequencies near the intrinsic frequency, the network dynamics are locked to the stimulation frequency (entrainment): this switch from intrinsic to driving frequency is mediated by an increase in forcing amplitude. The network was also found to exhibit subharmonic entrainment; that is, when the driving frequency is a harmonic of the intrinsic frequency, the network responds with a frequency that is subharmonic to the driving frequency and closest to the intrinsic frequency. This phenomenon has been observed previously in human EEG data when subjects were stimulated with different frequencies of flickering light and so-called steady-state visual evoked potentials (SSVEPs) were recorded (Herrmann, 2001). For example, at a stimulation frequency of 80 Hz, the visual cortex responded with a 10 Hz SSVEP oscillation; that is, a 1/8 subharmonic to the driving frequency. Therefore, and despite the simplicity of our model, our findings are consistent with recent tACS experiments showing that sinusoidal electric current applied to the human scalp resulted in entrainment of ongoing alpha oscillations (Helfrich et al., 2014). Our results also show good accordance with findings in animal models in which multiple driving frequencies were found to recruit the underlying Arnold entrainment tongues (Ali et al., 2013).
In contrast, our results indicate that the response dynamics to high-frequency stimulation (i.e., >50 Hz) builds on a different mechanism: network synchronous oscillations were found to be accelerated with respect to baseline through a nonlinear shift of the intrinsic frequency. Scaled by the stimulation intensity, the effect was characteristic of responses with amplified recurrent interactions, supporting a synergetic stability transition known to affect driven nonlinear networks generically (Lefebvre et al., 2013, 2015). Our result thus extend previous theoretical work (Lefebvre et al., 2015) and demonstrate that oscillatory nonlinear shifts can also be observed in sparse networks of spiking neurons using high-frequency stimuli. These computational predictions are also consistent with experimental findings showing such increases in peak alpha frequency during stimulation to flicker flashes of high frequency (Herrmann, 2001; Roberts and Robinson, 2012; cf. Fig. 1A). In addition, our results support the idea that high-frequency, sustained stimulation of cortical networks can shape the intrinsic alpha peak frequency and power via an alternative mechanism that could constitute a substitute and/or complementary strategy to engage cyclic neural dynamics.
At first glance, one might argue that our results do not compare to numerous experiments investigating direct cortical or transcranial brain stimulation because we used impulse trains rather than sinusoidal currents. It has been shown, however, that stimulating the visual cortex with either sinusoidally or rectangularly modulated light both resulted in sinusoidal SSVEPs (Dreyer and Herrmann, 2015). It is therefore plausible to assume that our findings can be used to better understand the human and animal findings on electric stimulation of brain tissue. In addition, our results replicate and extend the physiological findings of a recent SSVEP study that found a 1:1 Arnold tongue in the vicinity of the intrinsic frequency of ∼10 Hz in the human EEG (Notbohm et al., 2016). However, the number of parameter variations is limited in human and animal experiments, so network simulations that scan a broader range of stimulation frequencies and intensities are necessary.
Many studies that have used tACS have reported after-effects of enhanced amplitudes of ongoing brain oscillations after the end of rhythmical stimulation. It has been argued that such after-effects are due to synaptic plasticity because they do not appear at short stimulation durations (Strüber et al., 2015). Our simulation study did not model synaptic plasticity and thus did not investigate poststimulation after-effects. Nonetheless, our model provided a rich array of potential responses that could recruit and modulate plastic processes differentially. Future work should shed light on the stimulation parameters that are required to achieve such after-effects, which are likely to be critical if transcranial brain stimulation is to be used in treatment of diseases with unbalanced brain oscillations.
Footnotes
This work was supported by the Natural Sciences and Engineering Research Council of Canada (J.L.), the Swiss National Science Foundation (Grant PZ00P1_148186 to S.I. and Grants 320030-149982 and the National Centre of Competence in Research project SYNAPSY, The Synaptic Bases of Mental Disease, Project 51AU40-125759 to M.M.M.), the European Research Council under the European Union's Seventh Framework Programme [(FP7/2007–2013)/ERC Grant 257253 to A.H.], and the German Research Foundation (Grant HE3353/8-1).
The authors declare no competing financial interests.
- Correspondence should be addressed to Prof. Jérémie Lefebvre, Experimental Psychology Laboratory, Department of Psychology, Krembil Discovery Tower, 60 Leonard Avenue, Toronto, Ontario M5T 2S8, Canada. jeremie.lefebvre{at}uhnresearch.ca
