# Effect of Stimulation Waveform on the Non-linear Entrainment of Cortical Alpha Oscillations

^{1}Deutscher Wetterdienst, Department FE12-Data Assimilation, Offenbach am Main, Germany^{2}Rotman Research Institute, Baycrest Health Sciences, Toronto, ON, Canada^{3}Krembil Research Institute, University Health Network, Toronto, ON, Canada^{4}Experimental Psychology Lab, Department of Psychology, Cluster of Excellence “Hearing4all”, European Medical, School, Carl von Ossietzky University, Oldenburg, Germany^{5}Department of Mathematics and Institute for Biomaterials and Biomedical Engineering, University of Toronto, Toronto, ON, Canada

In the past decade, there has been a surge of interest in using patterned brain stimulation to manipulate cortical oscillations, in both experimental and clinical settings. But the relationship between stimulation waveform and its impact on ongoing oscillations remains poorly understood and severely restrains the development of new paradigms. To address some aspects of this intricate problem, we combine computational and mathematical approaches, providing new insights into the influence of waveform of both low and high-frequency stimuli on synchronous neural activity. Using a cellular-based cortical microcircuit network model, we performed numerical simulations to test the influence of different waveforms on ongoing alpha oscillations, and derived a mean-field description of stimulation-driven dynamics to better understand the observed responses. Our analysis shows that high-frequency periodic stimulation translates into an effective transformation of the neurons' response function, leading to waveform-dependent changes in oscillatory dynamics and resting state activity. Moreover, we found that randomly fluctuating stimulation linearizes the neuron response function while constant input moves its activation threshold. Taken together, our findings establish a new theoretical framework in which stimulation waveforms impact neural systems at the population-scale through non-linear interactions.

## Introduction

Oscillatory brain activity results from the collective and synchronous discharge of large populations of neurons, and is thought to play an important role in homeostasis, neural communication and information processing (Singer and Gray, 1995; Engel and Singer, 2001; Varela et al., 2001; Lakatos et al., 2008). In humans, such oscillations have been shown to be important for cognitive functions, and disturbed brain oscillations can result in cognitive deficits or neurological and psychiatric diseases (Uhlhaas and Singer, 2006). Many years of correlational analysis have shown that parameters of brain oscillations correlate with human perception, attention, memory, and behavior (Engel et al., 2001; Buzsáki and Draguhn, 2004; Hipp et al., 2011). Recent studies using TMS and tACS to modulate brain oscillations revealed a *causal* role of brain oscillations for such cognitive functions (e.g., Helfrich et al., 2014; Cecere et al., 2015; Dreyer and Herrmann, 2015). Importantly, all parameters of brain oscillations (amplitude, frequency, and phase) have been related to certain aspects of cognitive functions. While it has been repeatedly shown that e.g., tACS can up-regulate the amplitude of brain oscillations (Thut and Miniussi, 2009; Helfrich et al., 2014), less is known about down-regulating their amplitude or frequency, and how this depends on stimulation waveform. Intuitively, repetitive trains of negative and positive current pulses should have opposite effects on the frequency of brain oscillations. Do positive and negative pulse trains simply mirror each other with respect to the entrainment of alpha oscillations? Are pulses equivalent to sinusoids? Can noise-induced-like effects be triggered by deterministic signals? Answering these key questions would significantly improve our understanding of the role played by stimulation pattern on oscillatory brain dynamics, and catalyze the development of new clinical stimulation paradigms meant to engage neural populations and cortical oscillations.

This study sets out to answer some of these questions by harnessing computational and mathematical techniques and study the effect of stimulation waveform on cortical alpha oscillations. Alpha oscillations have been implicated in a wide variety of physiological and cognitive functions (Başar, 2012; Mierau et al., 2017), and have repeatedly been targeted using non-invasive stimulation in investigations aimed at obtaining a better understanding of the functional properties of cortical circuits (Fröhlich and McCormick, 2010; Cecere et al., 2015; Romei et al., 2016). Alpha oscillations have been shown to be maintained by large scale processes (Hindriks et al., 2014) supported by delayed network interactions (Cabral et al., 2014), and to build on slower and more global inhibitory processes (Klimesch et al., 2007; Womelsdorf et al., 2014). As such, to understand the effect of stimulation on these collective oscillations, a population-scale approach—in which networks of neurons are considered as opposed to individual cells—is necessary.

To provide new insight into the effects of brain stimulation on neural populations, we use two computational models in parallel and explore the impact of stimulation waveform and polarity on alpha oscillations. The first model, which we study numerically, is a cortical microcircuit network model which has been used before by the authors to investigate alpha resonance and entrainment in the cortex (Herrmann et al., 2016). The second model is a reduced neural oscillator model (Lefebvre et al., 2015; Hutt et al., 2016), which we derive from the cortical microcircuit model and analyze to understand the relationship between stimulation waveform and peak oscillation frequency. We combine insights provided by these two models to better understand, from a population-scale perspective, how stimulation waveforms can be tuned to either accelerate or slow down cortical alpha activity. We examine both near-resonant stimulation frequencies, where phase locking with the stimulation waveform can be observed, and higher frequencies, where we see the occurrence of non-linear entrainment. Through this approach, we develop a framework in which the effects of high-frequency stimuli with various waveform shapes on neural oscillations can be characterized by analyzing the associated transformation of the neurons' input/output (i.e., response) function. Recent experimental and computational studies have shown that the shape of the neural response function is altered in the presence of direct cortical stimulation, through a combination of somatic and synaptic effects (Lafon et al., 2017). To validate these results, and see how they are impacted by stimulation waveform, we systematically analyze neural population dynamics in the presence of repetitive pulse trains of positive and negative polarities, as well as sinusoidal drive and Gaussian white noise. We compare each case by deriving the associated mean-field dynamics, using a formalism that directly incorporates the effects of stimulation into the model equations. We then explore the influence of stimulation frequency and amplitude on network oscillations. In addition, our analysis suggests that, from a population perspective, the high-dimensionality of stimulation waveform parameter space can be significantly reduced by observing that seemingly distinct waveforms may possess equivalent entrainment properties. Our results further provide new perspectives on the waveform-specific interaction between stimulation and non-linear feedback in neural networks.

## Materials and Methods

### Cortical Microcircuit Model

To study the influence of different stimulation waveforms on oscillatory dynamics in cortical microcircuits, we here consider a model of interacting cortical populations and investigate changes in limit cycle solutions when subjected to stimulation. This model has been thoroughly discussed and analyzed in previous work (Herrmann et al., 2016), and thus we present it here briefly only.

This cortical network consists of spatially extended excitatory (*e*) and inhibitory (*i*) populations, whose activities are governed by the dynamics and interactions of neuronal ensembles. These ensembles, or sub-networks, include recurrently coupled neurons subjected to excitatory and inhibitory synaptic input, respectively. The ensemble spiking activity of each patch is modeled by the non-homogeneous Poisson processes

where ${{X}}_{{n}}^{{j}}{(}{t}{)}{=}{\displaystyle {\sum}_{{{t}}_{{l}}}{{\delta}}_{{n}}^{{j}}}{(}{t}{-}{{t}}_{{l}}{)}$ is the ensemble spike train of the *jth* patch and *n* = *e, i* indicate excitatory and inhibitory populations, respectively. The firing rate function*f*[*u*], also called the response function, sets the relationship between input potentials and output firing rates (Hutt and Buhry, 2014; Lefebvre et al., 2015; Herrmann et al., 2016). It exhibits a sigmoidal shape given by $f\left[{u}_{n}^{j}\right]={\left(1+exp\left[-\beta \left({u}_{n}^{j}-h\right)\right]\right)}^{-1}$ i.e., the firing rate probability approaches *f* = 1 for large membrane potentials, where the gain is β > 0 and the firing rate threshold is *h*. This defined, the model combines both the spiking of single cells as well as a dependence on the firing rate of the whole population. Such a hybrid cortical model thus combines both spiking and rate driven dynamics. The excitatory and inhibitory potentials ${u}_{e}^{j}\left(t\right)$ and ${u}_{i}^{j}\left(t\right)$ represent ensemble-averaged potentials proportional to averaged dendritic currents. They obey the set of non-linear stochastic equations

with *n* = *e, i* and the temporal rate constants α_{n}. The linear operator *L*[*U*] = *kU* represents membrane leaks. All cells in the network are driven by an external global stimulation *S*(*t*) exhibiting various waveforms. The cross-population recurrent inputs ${G}_{nm}^{j}\left(t\right)$ are defined by

where $PS{P}_{m}^{k}\left(t\right)$ refers to mean post-synaptic potential of patch *k* in population *m* at time *t*. Interactions between subnetworks *k* and *l* are subjected to intracortical propagation delaysτ^{jk} = |*x*(*k*)−*x*(*i*)| *v*^{−1}, with *v* being the axonal conduction velocity, set here to *v* = 0.13 m/s (Hutt et al., 2003). They are computed by convolving the time-delayed ensemble spike trains with exponential synapses of the form

with synaptic time constant *a*_{m}.

Excitatory and inhibitory populations are subjected to endogenous sources of noise ${\xi}_{n}^{j}\left(t\right)$, assumed to follow spatially and temporally independent Gaussian white noise profiles with fixed variance *D*_{n}. Synaptic weights within (${W}_{ee}^{jk}\left(c\right),{W}_{ii}^{jk}\left(c\right)$) and between (${W}_{ei}^{jk}\left(c\right),{W}_{ie}^{jk}\left(c\right)$) excitatory and inhibitory populations exhibit sparse exponential profiles (Hellwig, 2000) with connection probability *c*, that is

Neuron ensembles in the network are distributed randomly within a one-dimensional spatial domain Ω. The constants $\text{}{\sigma}_{n,m}^{2}={\sigma}_{e}^{2},\text{}{\sigma}_{i}^{2}$ correspond to the range of the excitatory and/or inhibitory interactions, *x*(*k*) refers to the spatial location of neurons in patch *k* and the connection probability is *c* = 0.6 i.e., 40% of the synaptic weights were randomly set to zero. The spatially-averaged neuroelectric network activity is a weighted sum over potentials of the excitatory and inhibitory population

where ${\varphi}_{e,i}^{k}$ are real positive coefficients. Here we assume that the network fine scale structure is unknown, and thus consider random weights i.e., ${\varphi}_{e,i}^{k}=$[0,1] (Herrmann et al., 2016). We did this to take into account various sources of observational variability that we do not model explicitly. However, specific choices of coefficient distributions can be made to increase the similarity of the neuroelectric output to signals such as LFPs and EEG (e.g., see Lindén et al., 2010). Model parameters are given in Table 1.

### Spectral Analysis

Spectral analysis was performed using a fast Fourier transform routine using freely available C++ scripts (Press et al., 2007). The power spectrum for each simulation condition is an average over five independent trials, each computed as the magnitude of the Fourier transform (with rectangular time window) of a time series of 4,000 ms duration. The long duration of the time series ensures negligible spectral leakage effects.

### Reduced Neural Oscillator Model

To better understand the mechanism involved in shaping oscillations in the cortical microcircuit model, we use a scalar and reduced non-linear network as a prototype to rigorously analyze the role of delayed and non-linear interactions in shaping emergent oscillations, and specifically how those are impacted by stimulation waveform. This simplified model sacrifices many physiological details in comparison to the cortical microcircuit model but preserves key components underlying the rhythmic activity seen in the cortical microcircuit model while remaining analytically tractable. Our goal here is to obtain a qualitative assessment of the different phenomena observed in our results.

Oscillations in the cortical microcircuit model arise due to delayed recurrent inhibition conveyed by inhibitory synapses. In this regime, inhibitory interactions dominate the dynamics, and the cortical microcircuit model can be significantly simplified, preserving the key components responsible of the oscillations. Specifically, we focus on parameters that result in an inhibition driven regime in which

Consequently, the dynamics of the cortical model obeys in good approximation

This approximation renders independent the dynamics of the inhibitory population from the activity of the excitatory cells. The excitatory membrane potential is thus, on average, driven by the activity of the inhibitory population, such that one may fully characterize the activity of the network by considering inhibitory ensemble dynamics. Assuming that the firing rate is high and that σ_{i} is small enough, i.e., broad spatial connectivity, we can write

where $\stackrel{\u0304}{\tau}=\underset{0}{\stackrel{\infty}{\int}}\tau g\left(\tau \right)d\tau \approx 25$ ms is the mean propagation delay and $g\left(\tau \right)=\frac{2\text{}{c}^{2}}{{\Omega}^{2}}\left(\frac{\Omega}{c}-\tau \right)|\text{}0\le \tau \le \frac{\Omega}{c}$. is the distribution of delays in our model. Taken together, now we may fully describe the response of the cortical populations to stimulation in the inhibitory-driven regime by the following scalar equation,

where ${\alpha}^{-1}\equiv {\alpha}_{i}^{-1}$ and ${U}^{j}\left(t\right)\equiv {u}_{i}^{j}\left(t\right)$ was introduced to distinguish the more detailed cortical microcircuit model and the reduced model in the subsequent calculations. The mean synaptic connectivity ${w}_{ii}^{o}=g{N}^{-1}$ has also been introduced to indicate the average evaluated over all possible pairs of inhibitory neurons. According to the derivations above the mean synaptic action is inhibitory with *g* < 0. The non-linear response function above remains the same with *f*[*u*] = (1+exp[−β(*u*−*h*)])^{−1}. This kind of approximation has been used frequently in the literature (e.g., Curtu and Ermentrout, 2004) to express the dynamics of excitatory and inhibitory networks from the perspective of a particular cellular species. For the rest of the analysis, we assume that *L*[*U*] = −*U*.

### Mean Field Dynamics in Presence of Stimulation

A common approach when trying to understand the essential dynamical characteristics of an otherwise high-dimensional system is to derive mean-field representations. What is different here is that we apply the mean field reduction by including stimulation in the calculations. As such, let us further assume that limit cycle solutions occur in a mean-driven regime in which the local dynamics can be seen as small independent fluctuations around a slowly varying mean Ū i.e.,

where Ū is given by

and < >_{N} is an average performed over the *N* units of the network. As an ansatz, local fluctuations *V*^{j} from the mean obey the zero mean processes

where we have used the fact that *L*[*U*] = −*U* and where

for *T* sufficiently small. Then taking the mean over *N* neurons in Equation (10) above yields the mean dynamics of the network in presence of stimulation

Now the network dynamics are governed by the effective neuron response function (Hutt et al., 2016)

where ρ(*V*) is the probability density function of the solution of Equation (13). According to this framework, the effect of a dynamic stimulus *S*(*t*) with stationary statistics can be characterized by looking at the probability density function ρ(*V*) associated to the linear and zero-mean processes whose dynamics obey Equation (13) and its convolution with the response function of the network as per Equation (16). We note that the mean field equation must be interpreted.

## Results

### Response of Cortical Neurons to Stimulation

To better understand the role of stimulation waveforms on the entrainment of network oscillations, we integrated numerically Equation (2) for different functional forms of the input term *S*(*t*): positive pulses, negative pulses, sinusoidal stimulation, and Gaussian white noise. In each case, all model parameters, except those related to the stimuli, were kept constant. Representative network responses to different waveforms are plotted in Figure 1. As seen in Figure 1A, without any stimulation, and for the set of parameters chosen, the network stabilizes into alpha-like synchronous activity. Spiking of the neurons is locked to these emergent global oscillations, resulting in a clear peak frequency of 10 Hz. When a pulse train stimulus with positive polarity was applied continuously at a rate of 50 Hz, network oscillations were found to accelerate with respect to baseline, stabilizing at a frequency of about 12 Hz. As shown in Figure 1B, the power of the associated oscillations and spike coherence were also both increased. We here recover the results of Herrmann et al. (2016), in which high-frequency (positive) pulse trains trigger non-linear acceleration of endogenous oscillations by changing the natural frequency of the solution of Equation (2). In Figure 1C however, when the network is stimulated continuously with a pulse train of negative polarity at a rate of 50 Hz, the opposite occurs: endogenous oscillations are slowed-down with respect to baseline. This novel effect, in contrast to the positive pulses, was not predicted by previous theoretical work. When sinusoidal stimulation was applied as shown in Figure 1D, endogenous oscillations were found to be entrained in a similar fashion as with positive pulses (cf. Figure 1A). Uncorrelated Gaussian white noise was found to have an analogous yet more pronounced effect (Figure 1E). Taken together, these results indicate that different stimulation waveforms have variable impact on the entrainment of endogenous network oscillations.

**Figure 1**. Variability of the cortical microcircuit network responses to stimulation with diverse waveforms. Stimulation has different impact on network resting state oscillations depending on the waveform applied. **(A)** Baseline oscillatory activity in absence of stimulation. The network neuroelectric output is plotted (top), as well as the spiking activity of the excitatory neurons (middle) and the power spectral density of the network output (bottom). The network displays stable oscillations of 10 Hz and ensemble spiking is phase locked to these intrinsic oscillations. **(B)** Positive pulses with a frequency of 50 Hz are found to accelerate endogenous oscillatory activity from 10 Hz **(A)** up to 12 Hz, entraining the neurons into faster cycles. Here, *S* = 2.5. Power expressed at the peak frequency also decreases with respect to baseline. **(C)** Negative pulses at 50 Hz, in contrast, decelerate endogenous oscillations from 10 Hz to about 8 Hz, and significantly suppress the power of endogenous oscillations. Here *S* = 2.5. **(D)** A sinusoidal waveform delivered at 50 Hz is also found to accelerate endogenous oscillations. The intensity of the sinusoidal signal, by virtue of being a continuous (i.e., not a discontinuous) was set to a smaller amplitude. Here *S* = 0.5. **(E)** For reference, Gaussian white noise (GWN) is also applied, and the 10 Hz oscillations are significantly accelerated from 10 Hz **(A)** up to 13 Hz. Here the noise intensity was increased from *D* = 0.0001 **(A–D)** to *D* = 0.01.

To understand how these results depend on stimulation settings and waveforms, we measured the response of the network while stimulation parameters were changed. For pulse trains and sinusoidal inputs, frequencies were systematically varied between 0 and 100 Hz with fixed amplitude. In the case of Gaussian white noise, the intensity of the noise was gradually increased between 0 and 0.01. Results are shown in Figure 2. Figure 2A shows that for positive pulses, increasing the stimulation frequency gradually shifted the peak response frequency from 10 to 12 Hz (Herrmann et al., 2016). The results were quite different with negative pulses: In Figure 2B, the peak frequency was slowed down and decreased in intensity (power decreases) until endogenous oscillations lost stability. Acceleration of endogenous oscillation was also observed in Figures 2C,D when sinusoidal input and Gaussian white noise were used, respectively. The different periodic stimulation cases are also depicted in Figure 3, where the range of stimulation frequencies is narrowed to a range about the endogenous frequency.

**Figure 2**. Diverse effects of stimulation waveform and frequency on the power spectrum in the cortical microcircuit network model. In each case, a waveform was chosen and used to stimulate the network. Stimulation frequency was increased while the power spectrum of the network responses was calculated. **(A)** Positive pulses, of frequencies ranging from 1 Hz up to 100 Hz are found to shape endogenous activity. The diagonal lines, representing the linear contribution of the stimulation waveform, delineate regions of entrainment. For slower frequencies, the network peak frequency shifts from its endogenous value to the stimulation's: neurons ensemble spiking is phased locked to the stimulation. For higher stimulation frequencies, endogenous oscillations accelerate, indicating non-linear entrainment. **(B)** Negative pulses do the opposite. Entrainment is weaker, and endogenous oscillations are gradually slowed down until they are suppressed as stimulation frequency is increased. **(C)** Sinusoidal inputs show a similar yet more pronounced effect as positive pulses. **(D)** By increasing the variance of the Gaussian white noise (GWN) input, the network's oscillations also accelerate. No entrainment can be seen as the stimulation possessed no dominant frequency. Parameters are as in Figure 1.

**Figure 3**. Impact of stimulation waveform on endogenous oscillations for near-resonant stimulation frequencies. This represents a close-up of the data plotted in Figure 2, over the interval 8–20 Hz. **(A)** Positive pulses with frequency ranging from 8 to 20 Hz shows that entrainment of network activity occurs for stimulation frequencies near the endogenous frequency of 10 Hz. The entrainment region is also non-symmetrical around that frequency: phase-locking is more pronounced to the right than to the left, indicating that the endogenous frequency accelerates as the stimulation frequency changes. **(B)** Negative pulses have a narrower entrainment region, and the network slowing down can readily be seen dominating the dynamics. **(C)** Sinusoidal stimulation is clearly more effective at entraining endogenous oscillations: the network dynamics is fully phase-locked to the stimulation frequency over this interval. Gaussian white noise (GWN) is here omitted because no entrainment occurs (see Figure 2). In all panels, the stimulation amplitude is the same. Parameters are as in Figure 2.

For stimulation frequencies close to but larger than the endogenous frequency, positive pulses and sinusoidal stimulation entrain the endogenous rhythm. Conversely, negative pulses entrain the endogenous rhythm for a more narrow range of frequencies.

### Impact of Stimulus Waveform on Alpha Oscillations: Theoretical Insights

To understand the mechanism behind the numerical observations made with the cortical microcircuit model (Figures 1, 2), we developed a reduced non-linear network model based on a mean-field approximation that preserved the mean features of the initial model, but remained analytically tractable (see section Materials and Methods).

As a first step to understand how the stimulation waveform affects endogenous oscillations, we applied the reduced neural oscillator model to characterize the effect of different stimulation waveforms on the response function of the network in Equation (15). These computations show that different stimulation patterns—leading to different statistics of the fluctuations around the activity mean—shape the effective response function in a plurality of waveform-dependent ways. The cases analyzed below are sequentially illustrated in Figure 4.

**Figure 4**. Effect of stimulation polarity and fluctuation distribution on the effective response function of the reduced model. Stimulation-driven fluctuations around the mean change the effective response function of the system in a waveform-dependent way. **(A)** In absence of stimulation, or in the case of a constant DC shift, the system's response function remains unchanged (top): the probability density function of the fluctuations is a delta function centered at 0. **(B)** For positive pulses, the network response function is linearized but primarily over the positive U axis. This is due to the asymmetrical shape of the probability density function. Note the change in effective threshold, due to the non-zero mean of the stimulation. Here *S* = 1, *r* = 0.2. **(C)** Analogous effect in presence of negative pulses, where the effective response function of the system is now linearized toward the negative U axis. Here *S* = −1, *r* = 0.2. **(D)** The sinusoidal input shapes the response function in a hybrid but symmetrical way, as its mean is centered at *U* = 0. Here *S* = 0.1, ω_{s} = 12 Hz. **(E)** Gaussian white noise (GWN) fully smoothes the system's response function, where the probability density function of the fluctuations are normally distributed and also centered around 0. Here *D* = 0.0001. All other parameters are as in Table 2.

**1. Pulse Train Stimulation**

Let us first consider the stimulus waveform

where δ(0) = 1 and zero otherwise. This pulse train has a rate *r* and an intensity of*S*. Using the mean-field formalism we have detailed earlier, fluctuations about the mean network activity ū obey

with μ_{S} = *Sr*. To obtain the effective neuron response function, one has to convolve the firing rate function *f* with the normalized probability density ρ(*V*) for the process (18) according to Equation (16)

In the limit of high gain, *f*[*U*] ≈ *H*[*U*−*h*], then the effective response function becomes

Then the network mean activity ū obeys the mean-field dynamics

Note that in these calculations, we have made no assumptions on the value of *S* and as such the result above holds for both positive and negative pulse trains. These two cases (*S* > 0, *S* < 0) are illustrated in Figures 4B,C.

**2. Sinusoidal Stimulation**

Let us now consider the periodic stimulation

with angle frequency ω_{s}. Fluctuations around the mean obey

since μ_{S} = 0. This case was studied in detail in Hutt et al. (2016). One can show that the associated normalized probability density function ρ(*V*) reads (Baker, 2006)

and consequently the effective non-linearity is given by

As illustrated in Figure 4D, the effective response function of the system is primarily linearized locally near the inflection point (i.e., *h* = 0).

**3. Gaussian White Noise Stimulation**

Next, we study Lefebvre and Hutt (2013):

where ξ_{j} are Gaussian white noise processes such that < ξ_{j} ξ _{k}> = δ_{jk} i.e., all neurons in the network experience independent stochastic input of intensity*D*. This case was also studied in detail in Hutt et al. (2016). Fluctuations around the mean obey the stochastic Langevin equations

where μ_{S} = 0. The associated probability density function ρ(*V*) is a symmetric Gaussian

Convolving with the threshold non-linearity yields

This equation states that additive noise linearizes the effective neuron response function as *D* is increased, cf. Figure 4E.

**4. Continuous Stimulation**

At last, let us consider the simple tonic stimulus

In this particular case, the input has zero variance and the effective response function is thus unchanged

The effect on the dynamics can be understood by introducing the change of variable

which leads to the mean–field dynamics

with $\stackrel{~}{F}\left[U\right]=H\left[U-\stackrel{~}{h}\right]$ and $\stackrel{~}{h}=h-S$. As such, the impact of continuous stimulation is analogous to a change in the response function threshold. This is in good agreement with the reported change in input/output response functions measured experimentally using tDCS (Lafon et al., 2017).

### Impact on Alpha Peak Frequency

As the derivations above have pointed out, stimulation statistics are reflected by changes in the effective response function, leading to mean-field equations with variable non-linear structures. Moreover, linear stability of the equilibria will also depend on stimulation statistics, which will be reflected on the features of oscillatory solutions. Regardless of stimuli waveform, limit cycle solutions are deployed around the implicitly defined equilibrium

The linearized dynamics around the steady state ū_{o} is

where $R\left[{\u016a}_{o}\right]=g\text{}{F}^{{\text{}}^{\prime}}\left[{\u016a}_{o}\right]$ and *R*[Ū_{o}] < 0 for the cases in Figure 4. We highlight here that the linear gain *R*[Ū_{o}] depends explicitly on the stimulation waveform through the convolved statistics in the effective function *F* and implicitly through the stimulus corrected equilibrium state Ū_{o}. To determine the frequency of emergent alpha oscillations, the iterative Galerkin method (He, 2005; Liu, 2005) is helpful to find a frequency ω minimizing the measure

Here we may consider small input stimuli i.e., *S*≪1and gain in the linear case

Using the ansatz Ū(*t*) = *A*cos(ω*t*)+Ū_{o}, one obtains for the first iteration

Setting *J*_{lin} = 0 and solving for ω yields an approximation of the linear frequency as a function of *R*[ū_{o}] and delay τ i.e.

Equation (36) approximates well the dependence of the network peak frequency on stimulation statistics whenever the stimulation amplitude remains small and its frequency high. Figure 5 illustrates how the network endogenous frequency depends on the linear gain for the delay considered in our model (i.e., $\stackrel{\u0304}{\tau}=25$ms). Although *R*[Ū_{o}] cannot always be computed analytically due to the implicit condition for the equilibrium (34), some specific cases such as Gaussian white noise for instance, remain tractable and accurate (Hutt et al., 2016), and can otherwise be computed numerically.

**Figure 5**. Peak frequency of the reduced model as a function of the linear gain. For a fixed delay, changes in the linear gain due to the stimulation will mediate the non-linear entrainment of the endogenous oscillations. Increases in |*R*| will cause a slowing down of endogenous oscillations, e.g., negative pulses, while decreases in |*R*| will do the opposite and increase the network peak frequency.

According to the local approximations above, effects of stimulation waveform on equilibrium states and oscillations can be characterized by local, stimulus-induced changes in the linearized gain, evaluated at the fixed point Ū_{o}. Equation (39) above states that ω is inversely proportional to |*R*|: increases (resp. decreases) in the linear gain *R* in Equation (35), i.e., decreases (resp. increases) of the slope of *F*[Ū_{o}], translates into an acceleration (resp. deceleration) of the network peak frequency. For the stimulation types studied and shown in Figure 4, this implies that network oscillations slow down whenever the system becomes locally non-linear, i.e., the transfer function becomes steeper, and accelerate when the network is pushed toward the linear regime (Hutt et al., 2016), i.e., the transfer function becomes more flat. Figure 4 (bottom panels) shows the non-linear response function and we observe that the slopes ${F}^{{\text{}}^{\prime}}\left[{\u016a}_{o}\right]$ of symmetric probability densities ρ are symmetric with respect to the threshold *h*, set here to 0, whereas the slopes of non-symmetric probability densities are non-symmetric. Since ${F}^{{\text{}}^{\prime}}\left[{\u016a}_{o}\right]~R$ for symmetric probability densities, the effect of external stimulation depends on the distance of the equilibrium state to the threshold only in contrast to non-symmetric probability densities. This highlights the importance of the shape of the probability density function and the equilibrium state.

According to this framework, mathematically the sensitivity of network oscillations to the stimulation waveform depends fully on how the probability density ρ(*V*) interacts with the network response function *f*. An important implication of this result is that different stimulus waveforms that possess the same statistical properties, i.e., the same probability density function ρ(*V*), will affect oscillations in the same way since the effective response function will possess an identical non-linear structure. This significantly reduces the dimensionality of the stimulus waveform parameter space, in the context of the optimization of non-linear acceleration/deceleration.

Using the approach outlined above, we computed the peak frequency of oscillations for various stimulation waveforms and compared them to the values computed numerically in the reduced network model. Results are presented in Figure 6. For all waveforms considered, the convolution approach (Equation 39) captured well the major effect of high-frequency stimuli on non-linear oscillations for small values of *S*, linking statistics of the input waveforms to the system's response. The mean field analysis confirms that positive (resp. negative) pulse trains accelerate (resp. slow down) ongoing oscillations, while both sinusoidal and Gaussian white noise inputs increased the endogenous frequency of the network.

**Figure 6**. Waveform-dependent non-linear entrainment of endogenous oscillations in the reduced model. The convolution approach allows us to compute the dependence of the system's peak frequency on the stimulation parameters. **(A)** Peak endogenous frequency as a function of the amplitude of pulses, either positive or negative (top). The theoretical calculations (red line), based on local linear analysis, are accurate for small pulse amplitudes, and diverge away from the simulated dynamics (black line) from there. The effect of pulses on the network peak frequency is monotonic and asymmetrical. The effect of the pulses' amplitude on the steady state ū_{o} is also plotted (bottom). **(B)** Expectedly, sinusoidal waveforms have a symmetrical effect whether the amplitude is positive or negative, and are found to accelerate the system's peak frequency in both cases. **(C)** Gaussian White noise (GWN) also accelerates endogenous oscillations, in a more pronounced way. **(D)** For comparison, a constant DC-shift input has an analogous effect as positive pulses, but the effect on the system's steady state is different.

The results above hold for high-frequency (or stochastic) stimulation; what happens if this condition is relaxed? We investigated this question numerically and results are plotted in Figure 7. For slower frequencies, the network dynamics were found to be dominated by (super- and sub-) harmonic entrainment: endogenous oscillations maintain stable phase relationships with the stimuli, sequentially jumping from one harmonic to the next, as seen from the Arnold tongue patterns portrayed in both Figures 7A,B. There, the waveform-dependent non-linear interactions between stimulation and network oscillations can be observed through the presence of tilted Arnold tongues, due to non-linear shifts in the natural frequency of the system (Lefebvre et al., 2015; Herrmann et al., 2016; Hutt et al., 2016). Surprisingly, negative pulses were found to be more efficient at entraining ongoing oscillations. Indeed, Figures 7A,B show that Arnold tongues are wider for negative pule trains (*S* < 0) compared to positive ones (*S* > 0). Also, negative pulses were found to suppress network oscillations beyond specific values of stimulation amplitude and frequency. We note that this is fully analogous with the dynamics observed in Figure 2 where peak power is gradually suppressed by negative pulse trains whose frequency exceeds ~50 Hz for *S* = 2.5. Beyond this point, the network does not exhibit any internal resonances and the peak frequency is the same as the stimulation frequency.

**Figure 7**. Numerically computed Arnold Tongues of reduced model for different waveforms. Given the different impact of the stimulation on the network response function and oscillatory properties, different waveforms possess different entrainment properties (linear or not). **(A)** Peak endogenous frequency as a function of varying stimulation frequency and amplitude for pulse waveforms. The triangular regions represent Arnold tongues were the system's activity is phase locked to the stimulus drive. The resulting map is asymmetrical, where acceleration can be seen for positive pulses, and slowing down (and loss of stability) can be observed for negative pulses. Notice also that the Arnold tongues are bent to the right, a clear signature of a change in endogenous frequency due to non-linear entrainment. The dashed line denotes the parameter border below which network activity vanishes. **(B)** In contrast, sinusoidal inputs have fully symmetrical entrainment regions. These are analogously bent to the right, as expected.

## Discussion

Electromagnetic brain stimulation has become increasingly popular to support a wide variety of clinical interventions. It is routinely used in the treatment of various types of neuropsychiatric disorders such as treatment-resistant depression (Ferrucci et al., 2009), Parkinsonism (Hess, 2013), Schizophrenia (Hoy et al., 2016), and the number of potential applications is increasing rapidly. Both invasive (e.g., DBS) and non-invasive (e.g., TMS, TACS) paradigms are emerging as strong alternatives to pharmaceutical treatments and have further raised the fascinating prospect of entraining brain rhythms to engage neural circuits at a functional level (Fröhlich, 2015). However, the relationship between stimulation waveform and entrainment outcomes remains a difficult and open problem. How the different stimulation temporal patterns interact with the intrinsic non-linear structure of cortical circuits to shape synchronous neural dynamics is also poorly understood.

To answer some of these questions, we have combined computational and mathematical methods to reconcile the effect of stimulation at the mesoscopic neural ensemble and macroscopic population scales using mean-field techniques. Using a detailed cortical microcircuit model, we have numerically explored the effect of stimulation pulse trains, sinusoidal inputs and Gaussian white noise stimulation on resting state alpha oscillations. First, our simulations have confirmed that, as expected, distinct waveforms have different entrainment properties. Notably, while positive pulse trains were found to accelerate ongoing oscillations, negative pulses did the opposite. To understand the source of this novel finding, we have developed a framework in which the mesoscopic and waveform-dependent effects of stimulation on oscillatory activity can be characterized by a change in the neuron response function, mathematically speaking through a convolution with the probability density function associated with stimulation-induced fluctuations around the mean. Using this approach, it was possible to relate the statistics of stimuli to the acceleration/slowing down of non-linear oscillations using a reduced neural oscillator model that preserved the non-linear structure of the more detailed cortical microcircuit model. Taken together, these results show that the core differences in entrainment properties between various stimulation waveforms can be explained by population-scale changes in the effective response function of the network—an emerging perspective in line with recent findings in the literature (Kar et al., 2017; Lafon et al., 2017). The approach we have put forward is relevant especially in cases where the relative phase of stimulation with respect to ongoing oscillations is not known and where the stimulus and the ongoing activity coexist on different time scales. Brain stimulation, especially in the non-invasive case, engages populations of neurons with effects that go beyond single cells and extend whole circuits. As such, the use of mean-field techniques is warranted, providing appropriate description of the dynamics at these larger spatial scales. The non-linear structure of neuron response functions reflects the combined effect of multiple neurophysiological mechanisms such as ion channel dynamics (e.g., Hutt and Buhry, 2014), adaptation (e.g., Benda et al., 2007) as well as recurrent feedback (e.g., Sutherland et al., 2009). As such, our results suggest that distinct waveforms engage neural circuits through a plurality of neurophysiological mechanisms, yet generate equivalent outcomes as far as synchronous oscillations are concerned.

The analysis we have conducted has revealed that in regimes of high frequency stimulation, distinct waveforms may have similar impact on non-linear neural oscillations. Indeed, positive pulses, sinusoidal inputs and Gaussian white noise were all found to have similar influence on ongoing oscillations when the stimulation frequency (or intensity *D* for Gaussian white noise) is high, indicating that, within these regimes, it is possible to reduce the high-dimensionality of the stimulation waveform parameter space: multiple stimulation patterns have equivalent effect on synchronous oscillations.

Our analysis has shown that, in first approximation, the local features of the effective response function *F* near the fixed point Ū_{o} determine primarily the impact of stimulation waveforms on limit cycle solutions. This means that the stimulation efficacy is highly dependent on the state: changes to the response function are commensurate with the proximity of the fixed point, implying that the same waveform will have different impact on endogenous oscillations if the equilibrium location in phase space differs. This is in complete agreement with numerous recent findings showing the state dependence of entrainment efficacy (Neuling et al., 2013; Alagapan et al., 2016).

The results presented in this study are aimed at obtaining a better understanding of the properties of mesoscopic neural population activity in regimes where the dynamics are dominated by recurrent inhibition. In our model, this arises due to differential spatial profiles of excitatory vs. inhibitory connections. Specifically, excitatory connections dominate at smaller spatial scales, whereas inhibitory connections dominate at larger spatial scales. This is the 1-D analog of classic “Mexican Hat” profiles of lateral connectivity, as observed for example in visual cortex (Amari, 1977; Kang et al., 2003). Other neural population-scale mechanisms of rhythmogenesis that have been studied in the experimental and theoretical literature include intra-columnar circuit motifs (Wilson and Cowan, 1972; Jansen and Rit, 1995; Womelsdorf et al., 2014) and long-range thalamocortical loops (Lopes da Silva et al., 1974; Robinson et al., 2001). Distinguishing the differential contributions of these and other mechanisms to a given set of empirical observations is a challenging and active research topic. However, despite major differences in their exact neuroanatomical and neurophysiological components, many of these different circuit mechanisms can be understood as having in common the same core rhythmogenic phenomenon—namely delayed recurrent inhibition. As such, the mathematical techniques developed and insights obtained in the present study regarding changes in neural population response functions under different stimulation waveform regimes may, we suggest, also prove useful for characterizing analogous responses in other modeling and experimental contexts.

In order to fully characterize the effects of electromagnetic brain stimulation, it is necessary to consider both local effects at the stimulation site, and distal effects that result from propagation of stimulated response throughout the brain via long-range white matter pathways (Massimini et al., 2005). Our focus in this study has been on the first of these, with the neural population and mean field models presented describing oscillatory responses to various stimulation waveforms in an isolated patch of cortex. An important direction for future work shall be to extend this approach to investigate polysynaptic and large-scale network effects. Such networks, extending on larger spatial scales, are built of interconnected population patches, as considered here. The rapidly growing literature on stimulation and control problems in macro-connectomics has primarily used relatively simple descriptions of neural population dynamics such as planar oscillators (e.g., Spiegler et al., 2016), phase oscillators (e.g., Gollo et al., 2017), or linear integrators (e.g., Betzel et al., 2016), and typically focuses on spatial pattern formation and recovery or modulation of canonical fmri-derived resting state networks. A number of studies have examined power spectrum changes in M/EEG-relevant frequency ranges following focal stimulation (Spiegler et al., 2010; Cona et al., 2011, 2014; Fung and Robinson, 2014). Recent work by Spiegler and colleagues (Kunze et al., 2016) observed a sharpening and acceleration of peak oscillatory frequencies, as well as an increase in long-range synchrony, following simulated transcranial direct current stimulation in a model of coupled neural masses. Importantly, these authors found that heterogeneity in the level of neural population response function saturation—which arose due to the network topology and to the spatially diffuse effect of electrical brain stimulation—played a major role in shaping large-scale network activity. Likewise, developing a spatially mapped characterization of node-wise variation in response function modifications due to diffuse and inhomogeneous stimulation current distributions shall be important for translating the theoretical insights obtained in the present study to improved experimental paradigms and clinical treatments. Taken together and in conjunction with previous studies, our results reinforce the notion that stimulation effects and optimization has to be considered from multiple spatial scales. Indeed, previous work from the authors have demonstrated that delays, which become non-negligible as the spatial scale considered increases, play an important role not only on the immediate response of the stimulated networks (Hutt and Atay, 2007; Lefebvre et al., 2015; Hutt et al., 2016), but also on post-stimulation after effects (Alagapan et al., 2016), and thus play a key role in defining the entrainment properties of non-linear neural systems.

## Author Contributions

AH and JL conceived the principle idea of the work. JL performed the numerical analysis. AH, JG, and JL have worked out the analytical results. JL and CH structured the manuscript and all authors have written the manuscript.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This work has been supported by the Natural Sciences and Engineering Research Council of Canada (JL).

## References

Alagapan, S., Schmidt, S. L., Lefebvre, J., Hadar, E., Shin, H. W., and Fröhlich, F. (2016). Modulation of cortical oscillations by low-frequency direct cortical stimulation is state-dependent. *PLoS Biol.* 14:e1002424. doi: 10.1371/journal.pbio.1002424

Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. *Biol. Cybern.* 27, 77–87. doi: 10.1007/BF00337259

Baker, G. L. (2006). Probability, pendulums, and pedagogy. *Am. J. Phys.* 74, 482–489. doi: 10.1119/1.2186689

Başar, E. (2012). A review of alpha activity in integrative brain function: fundamental physiology, sensory coding, cognition and pathology. *Int. J. Psychophysiol.* 86, 1–24. doi: 10.1016/j.ijpsycho.2012.07.002

Benda, J., Bethge, M., Hennig, M., Pawelzik, K., and Herz, A. V. M. (2007). Spike-frequency adaptation: phenomenological model and experimental tests. *Neurocomputing* 38–40, 105–110. doi: 10.1016/S0925-2312(01)00545-8

Betzel, R., Gu, S., Medaglia, J., Pasqualetti, F., and Bassett, D. (2016). Optimally controlling the human connectome: the role of network topology. *Sci. Rep.* 6:30770. doi: 10.1038/srep30770

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. *Science* 304:1926. doi: 10.1126/science.1099745

Cabral, J., Luckhoo, H., Woolrich, M., Joensson, M., Mohseni, H., Baker, A., et al. (2014). Exploring mechanisms of spontaneous functional connectivity in MEG: how delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. *Neuroimage* 90, 423–435. doi: 10.1016/j.neuroimage.2013.11.047

Cecere, R., Rees, G., and Romei, V. (2015). Individual differences in alpha frequency drive crossmodal illusory perception. *Curr. Biol*. 19, 231–235. doi: 10.1016/j.cub.2014.11.034

Cona, F., Lacanna, M., and Ursino, M. (2014). A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. *J. Comput. Neurosci.* 37, 125–148. doi: 10.1007/s10827-013-0493-1

Cona, F., Zavaglia, M., Massimini, M., Rosanova, M., and Ursino, M. (2011). A neural mass model of interconnected regions simulates rhythm propagation observed via TMS-EEG. *Neuroimage* 57, 1045–1058. doi: 10.1016/j.neuroimage.2011.05.007

Curtu, R., and Ermentrout, B. (2004). Pattern formation in a network of excitatory and inhibitory cells with adaptation. *SIAM J. Appl. Dyn. Syst*. 3, 191–231. doi: 10.1137/030600503

Dreyer, A. M., and Herrmann, C. S. (2015). Frequency-modulated steady-state visual evoked potentials: a new stimulation method for brain-computer interfaces. *J Neurosci. Methods* 241, 1–9. doi: 10.1016/j.jneumeth.2014.12.004

Engel, A. K., and Singer, W. (2001). Temporal binding and the neural correlates of sensory awareness. *Trends Cogn. Sci.* 5, 16–25. doi: 10.1016/S1364-6613(00)01568-0

Engel, A. K., Fries, P., and Singer, W. (2001). Dynamic predictions: oscillations and synchrony in top-down processing. *Nat. Rev. Neurosci*. 2, 704–716. doi: 10.1038/35094565

Ferrucci, R., Bortolomasi, M., Vergari, M., Tadini, L., Salvoro, B., Giacopuzzi, M., et al. (2009). Transcranial direct current stimulation in severe, drug-resistant major depression. *J. Affect. Disord.* 118, 215–219. doi: 10.1016/j.jad.2009.02.015

Fröhlich, F. (2015). Experiments and models of cortical oscillations as a target for non-invasive brain stimulation. *Prog. Brain Res.* 222, 41–73. doi: 10.1016/bs.pbr.2015.07.025

Fröhlich, F., and McCormick, D. A. (2010). Endogenous electric fields may guide neocortical network activity. *Neuron* 67, 129–143. doi: 10.1016/j.neuron.2010.06.005

Fung, P. K., and Robinson, P. A. (2014). Neural field theory of synaptic metaplasticity with applications to theta burst stimulation. *J. Theor. Biol.* 340, 164–176. doi: 10.1016/j.jtbi.2013.09.021

Gollo, L., Roberts, J., and Cocchi, L. (2017). Mapping how local perturbations influence systems-level brain dynamics. *Neuroimage* 60, 97–112. doi: 10.1016/j.neuroimage.2017.01.057

He, J. H. (2005). Periodic solutions and bifurcations of delay-differential equations. *Phys. Lett.* A 347, 228–230. doi: 10.1016/j.physleta.2005.08.014

Helfrich, R. F., Schneider, T. R., Rach, S., Trautmann-Lengsfeld, S. A., Engel, A. K., and Herrmann, C. S. (2014). Entrainment of brain oscillations by transcranial alternating current stimulation. *Curr. Biol.* 24, 333–339. doi: 10.1016/j.cub.2013.12.041

Hellwig, B. (2000). A quantitative analysis of the local connectivity between pyramidal neurons in layers 2/3 of the rat visual cortex. *Biol. Cybern.* 82, 111–121. doi: 10.1007/PL00007964

Herrmann, C., Ionta, S., Hutt, A., Murray, M. M., and Lefebvre, J. (2016). Shaping intrinsic oscillatory neural activity with periodic stimulation. *J. Neurosci.* 36, 5328–5337. doi: 10.1523/JNEUROSCI.0236-16.2016

Hess, C. W. (2013). Modulation of cortical-subcortical networks in Parkinson's disease by applied field effects. *Front. Hum. Neurosci.* 7:565. doi: 10.3389/fnhum.2013.00565

Hindriks, R., van Putten, M. J. A., and Deco, G. (2014). Intra-cortical propagation of EEG alpha oscillations. *Neuroimage* 103, 444–453. doi: 10.1016/j.neuroimage.2014.08.027

Hipp, J. F., Engel, A. K., and Siegel, M. (2011). Oscillatory synchronization in large-scale cortical networks predicts perception. *Neuron* 69, 387–396. doi: 10.1016/j.neuron.2010.12.027

Hoy, K. E., Whitty, D., Bailey, N., and Fitzgerald, P. B. (2016). Preliminary investigation of the effects of γ-tACS on working memory in schizophrenia. *J Neural Transm.* 123, 1205–1212. doi: 10.1007/s00702-016-1554-1

Hutt, A., and Atay, F. M. (2007). Spontaneous and evoked activity in extended neural populations with gamma-distributed spatial interactions and transmission delay. *Chaos Sol. Fract*. 32, 547–560. doi: 10.1016/j.chaos.2005.10.091

Hutt, A., and Buhry, L. (2014). Study of GABAergic extra-synaptic tonic inhibition in single neurons and neural populations by traversing neural scales: application to propofol-induced anaesthesia. *J Comput Neurosci.* 37, 417–437. doi: 10.1007/s10827-014-0512-x

Hutt, A., Bestehorn, M., and Wennekers, T. (2003). Pattern formation in intracortical neuronal fields. *Netw. Comput. Neural Syst*. 14, 351–368. doi: 10.1088/0954-898X_14_2_310

Hutt, A., Mierau, A., and Lefebvre, J. (2016). Dynamic control of synchronous dynamics in networks of spiking neurons. *PLoS ONE* 11:e0161488. doi: 10.1371/journal.pone.0161488

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. *Biol. Cybern.* 73, 357–366. doi: 10.1007/BF00199471

Kang, K., Shelley, M., and Sompolinsky, H. (2003). Mexican hats and pinwheels in visual cortex. *Proc. Natl. Acad. Sci. U.S.A*. 100, 2848–2853. doi: 10.1073/pnas.0138051100

Kar, K., Duijnhouwer, J., and Krekelberg, B. (2017). Transcranial alternating current stimulation attenuates neuronal adaptation. *J. Neurosci.* 30, 2266–2216. doi: 10.1523/JNEUROSCI.2266-16.2016

Klimesch, W., Sauseng, P., and Hanslmayr, S. (2007). EEG alpha oscillations: the inhibition-timing hypothesis. *Brain Res. Rev.* 53, 63–88. doi: 10.1016/j.brainresrev.2006.06.003

Kunze, T., Hunold, A., Haueisen, J., Jirsa, V., and Spiegler, A. (2016). Transcranial direct current stimulation changes resting state functional connectivity: a large-scale brain network modeling study. *Neuroimage* 140, 174–187. doi: 10.1016/j.neuroimage.2016.02.015

Lafon, B., Rahman, A., Bikson, M., and Parra, L. C. (2017). Direct Current stimulation alters neuronal input/output function. *Brain Stimul.* 10, 36–45. doi: 10.1016/j.brs.2016.08.014

Lakatos, P., Karmos, G., Mehta, A. D., Ulbert, I., and Schroeder, C. E. (2008). Entrainment of neuronal oscillations as a mechanism of attentional selection. *Science* 320, 110–113. doi: 10.1126/science.1154735

Lefebvre, J., and Hutt, A. (2013). Additive noise quenches delay-induced oscillations. *Europhys. Lett*. 102:60003. doi: 10.1209/0295-5075/102/60003

Lefebvre, J., Hutt, A., Knebel, J. F., Whittingstall, K., and Murray, M. M. (2015). Stimulus statistics shape oscillations in nonlinear recurrent neural networks. *J. Neurosci*. 35, 2895–2903. doi: 10.1523/JNEUROSCI.3609-14.2015

Lindén, H., Pettersen, K. H., and Einevoll, G. T. (2010). Intrinsic dendritic filtering gives low-pass power spectra of local field potentials. *J. Comput. Neurosci.* 29, 423–444. doi: 10.1007/s10827-010-0245-4

Liu, H. M. (2005). Approximate period of nonlinear oscillators with discontinuities by modified Lindstedt-Poincare method. *Chaos Solit. Fract.* 23, 577–579. doi: 10.1016/j.chaos.2004.05.004

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of brain rhythmic activity. *Kybernetik* 15, 27–37. doi: 10.1007/BF00270757

Massimini, M., Ferrarelli, F., Huber, R., Esser, S. K., Singh, H., and Tononi, G. (2005). Breakdown of cortical effective connectivity during sleep. *Science* 309, 2228–2232. doi: 10.1126/science.1117256

Mierau, A., Klimesch, W., and Lefebvre, J. (2017). State-dependent alpha peak frequency shifts: experimental evidence, potential mechanisms and functional implications. *Neuroscience* 360, 146–154. doi: 10.1016/j.neuroscience.2017.07.037

Neuling, T., Rach, S., and Herrmann, C. S. (2013). Orchestrating neuronal networks: sustained after-effects of transcranial alternating current stimulation depend upon brain states. *Front. Human Neurosci*. 7:161. doi: 10.3389/fnhum.2013.00161

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2007). *Numerical Recipes: The Art of Scientific Computing, 3rd Edn*. Cambridge: Cambridge University Press.

Robinson, P. A., Rennie, C. J., Wright, J. J., Bahramali, H., Gordon, E., and Rowe, D. L. (2001). Prediction of electroencephalographic spectra from neurophysiology. *Phys. Rev. E* 63:021903. doi: 10.1103/PhysRevE.63.021903

Romei, V., Bauer, M., Brooks, J. L., Economides, M., Penny, W., Thut, G., et al. (2016). Causal evidence that intrinsic beta-frequency is relevant for enhanced signal propagation in the motor system as shown through rhythmic TMS. *Neuroimage* 125, 120–130. doi: 10.1016/j.neuroimage.2015.11.020

Singer, W., and Gray, C. M. (1995). Visual feature integration and the temporal correlation hypothesis. *Annu. Rev. Neurosci.* 18, 555–586. doi: 10.1146/annurev.ne.18.030195.003011

Spiegler, A., Hansen, E., Bernard, C., McIntosh, A., and Jirsa, V. (2016). Selective activation of resting-state networks following focal stimulation in a connectome-based network model of the human brain. *eNeuro* 3:ENEURO.0068–16.2016. doi: 10.1523/ENEURO.0068-16.2016

Spiegler, A., Kiebel, S. J., Atay, F. M., and Knösche, T. R. (2010). Bifurcation analysis of neural mass models: impact of extrinsic inputs and dendritic time constants. *Neuroimage* 52, 1041–1058. doi: 10.1016/j.neuroimage.2009.12.081

Sutherland, C., Doiron, B., and Longtin, A. (2009). Feedback-induced gain control in stochastic spiking networks. *Biol. Cybern.* 100, 475–489. doi: 10.1007/s00422-009-0298-5

Thut, G., and Miniussi, C. (2009). New insights into rhythmic brain activity from TMS–EEG studies. *Trends Cogn. Sci.* 13, 182–189. doi: 10.1016/j.tics.2009.01.004

Uhlhaas, P. J., and Singer, W. (2006). Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. *Neuron* 52, 155–168. doi: 10.1016/j.neuron.2006.09.020

Varela, F., Lachaux, J.-P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. *Nat. Rev. Neurosci.* 2, 229–239. doi: 10.1038/35067550

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. *Biophys. J.* 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

Keywords: stimulation waveform, synchrony, entrainment, neural dynamics, networks, oscillations

Citation: Hutt A, Griffiths JD, Herrmann CS and Lefebvre J (2018) Effect of Stimulation Waveform on the Non-linear Entrainment of Cortical Alpha Oscillations. *Front. Neurosci*. 12:376. doi: 10.3389/fnins.2018.00376

Received: 02 March 2018; Accepted: 16 May 2018;

Published: 26 June 2018.

Edited by:

Olivier David, Institut National de la Santé et de la Recherche Médicale (INSERM), FranceReviewed by:

Thomas R. Knösche, Max-Planck-Institut für Kognitions- und Neurowissenschaften, GermanyTamer Demiralp, Istanbul University, Turkey

Copyright © 2018 Hutt, Griffiths, Herrmann and Lefebvre. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jérémie Lefebvre, jeremie.lefebvre@uhnresearch.ca