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

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;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 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 noninvasive 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 populationscale 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 meanfield 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.

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 is the ensemble spike train of the jth patch and n = e, i indicate excitatory and inhibitory populations, respectively. The firing rate functionf [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 u 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 j e (t) and u j i (t) 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] = k U 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 j nm (t) are defined by where PSP k m (t) 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 ξ j n (t), assumed to follow spatially and temporally independent Gaussian white noise profiles with fixed variance D n . Synaptic weights within (W jk ee (c) , W jk ii (c)) and between (W jk ei (c) , W jk ie (c)) 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 σ 2 n,m = σ 2 e , σ 2 i 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 φ k e,i are real positive coefficients. Here we assume that the network fine scale structure is unknown, and thus consider random weights i.e., φ k e,i =[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 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, was introduced to distinguish the more detailed cortical microcircuit model and the reduced model in the subsequent calculations. The mean synaptic connectivity w o ii = gN −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 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 bȳ 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 ) 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.

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 nonlinear 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. 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.
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.

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 ofS. 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) ρ 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.

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 bỹ As illustrated in Figure 4D, the effective response function of the system is primarily linearized locally near the inflection point (i.e., h = 0).

Gaussian White Noise Stimulation
Next, we study Lefebvre and Hutt (2013):  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.  Table 2. where ξ j are Gaussian white noise processes such that < ξ j ξ k >= δ jk i.e., all neurons in the network experience independent stochastic input of intensityD. 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 (29) This equation states that additive noise linearizes the effective neuron response function as D is increased, cf. Figure 4E.

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 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 and R Ū o < 0 for the cases in 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 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.
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.,τ = 25ms). 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 , and can otherwise be computed numerically. 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 , 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 ′ [Ū o ] 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 nonsymmetric. Since F ′ Ū o ∼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 nonsymmetric 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 nonlinear 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.
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 waveformdependent 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.

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 stimulationinduced 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 meanfield 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 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 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.
(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 U 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., 2011Cona et al., , 2014Fung 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 longrange 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 largescale 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.