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

Abstract

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 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 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 and 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 are defined by

where 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⁻¹, 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 , assumed to follow spatially and temporally independent Gaussian white noise profiles with fixed variance D_n. Synaptic weights within () and between () 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 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 are real positive coefficients. Here we assume that the network fine scale structure is unknown, and thus consider random weights i.e., [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 ms is the mean propagation delay and . 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 and was introduced to distinguish the more detailed cortical microcircuit model and the reduced model in the subsequent calculations. The mean synaptic connectivity 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)])⁻¹. 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

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

Figure 3

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

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 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)

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

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 and . 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 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) = Acos(ω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., 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

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

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

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.

References

• 1

  AlagapanS.SchmidtS. L.LefebvreJ.HadarE.ShinH. W.FröhlichF. (2016). Modulation of cortical oscillations by low-frequency direct cortical stimulation is state-dependent. PLoS Biol.14:e1002424. 10.1371/journal.pbio.1002424

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BakerG. L. (2006). Probability, pendulums, and pedagogy. Am. J. Phys.74, 482–489. 10.1119/1.2186689

  • CrossRef
  • Google Scholar

• 4

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BendaJ.BethgeM.HennigM. K PawelzikK.HerzA. V. M. (2007). Spike-frequency adaptation: phenomenological model and experimental tests. Neurocomputing38–40, 105–110. 10.1016/S0925-2312(01)00545-8

  • CrossRef
  • Google Scholar

• 6

  BetzelR.GuS.MedagliaJ.PasqualettiF.BassettD. (2016). Optimally controlling the human connectome: the role of network topology. Sci. Rep.6:30770. 10.1038/srep30770

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BuzsákiG.DraguhnA. (2004). Neuronal oscillations in cortical networks. Science 304:1926. 10.1126/science.1099745

  • CrossRef
  • Google Scholar

• 8

  CabralJ.LuckhooH.WoolrichM.JoenssonM.MohseniH.BakerA.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. Neuroimage90, 423–435. 10.1016/j.neuroimage.2013.11.047

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  CecereR.ReesG.RomeiV. (2015). Individual differences in alpha frequency drive crossmodal illusory perception. Curr. Biol. 19, 231–235. 10.1016/j.cub.2014.11.034

  • CrossRef
  • Google Scholar

• 10

  ConaF.LacannaM.UrsinoM. (2014). A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. J. Comput. Neurosci.37, 125–148. 10.1007/s10827-013-0493-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  ConaF.ZavagliaM.MassiminiM.RosanovaM.UrsinoM. (2011). A neural mass model of interconnected regions simulates rhythm propagation observed via TMS-EEG. Neuroimage57, 1045–1058. 10.1016/j.neuroimage.2011.05.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  CurtuR.ErmentroutB. (2004). Pattern formation in a network of excitatory and inhibitory cells with adaptation. SIAM J. Appl. Dyn. Syst. 3, 191–231. 10.1137/030600503

  • CrossRef
  • Google Scholar

• 13

  DreyerA. M.HerrmannC. S. (2015). Frequency-modulated steady-state visual evoked potentials: a new stimulation method for brain-computer interfaces. J Neurosci. Methods241, 1–9. 10.1016/j.jneumeth.2014.12.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  EngelA. K.SingerW. (2001). Temporal binding and the neural correlates of sensory awareness. Trends Cogn. Sci.5, 16–25. 10.1016/S1364-6613(00)01568-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  EngelA. K.FriesP.SingerW. (2001). Dynamic predictions: oscillations and synchrony in top-down processing. Nat. Rev. Neurosci. 2, 704–716. 10.1038/35094565

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  FerrucciR.BortolomasiM.VergariM.TadiniL.SalvoroB.GiacopuzziM.et al. (2009). Transcranial direct current stimulation in severe, drug-resistant major depression. J. Affect. Disord.118, 215–219. 10.1016/j.jad.2009.02.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

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

  • CrossRef
  • Google Scholar

• 18

  FröhlichF.McCormickD. A. (2010). Endogenous electric fields may guide neocortical network activity. Neuron67, 129–143. 10.1016/j.neuron.2010.06.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  FungP. K.RobinsonP. A. (2014). Neural field theory of synaptic metaplasticity with applications to theta burst stimulation. J. Theor. Biol.340, 164–176. 10.1016/j.jtbi.2013.09.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  GolloL.RobertsJ.CocchiL. (2017). Mapping how local perturbations influence systems-level brain dynamics. Neuroimage60, 97–112. 10.1016/j.neuroimage.2017.01.057

  • CrossRef
  • Google Scholar

• 21

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

  • CrossRef
  • Google Scholar

• 22

  HelfrichR. F.SchneiderT. R.RachS.Trautmann-LengsfeldS. A.EngelA. K.HerrmannC. S. (2014). Entrainment of brain oscillations by transcranial alternating current stimulation. Curr. Biol.24, 333–339. 10.1016/j.cub.2013.12.041

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  HellwigB. (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. 10.1007/PL00007964

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  HerrmannC.IontaS.HuttA.MurrayM. M.LefebvreJ. (2016). Shaping intrinsic oscillatory neural activity with periodic stimulation. J. Neurosci.36, 5328–5337. 10.1523/JNEUROSCI.0236-16.2016

  • CrossRef
  • Google Scholar

• 25

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  HindriksR.van PuttenM. J. A.DecoG. (2014). Intra-cortical propagation of EEG alpha oscillations. Neuroimage103, 444–453. 10.1016/j.neuroimage.2014.08.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  HippJ. F.EngelA. K.SiegelM. (2011). Oscillatory synchronization in large-scale cortical networks predicts perception. Neuron69, 387–396. 10.1016/j.neuron.2010.12.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  HoyK. E.WhittyD.BaileyN.FitzgeraldP. B. (2016). Preliminary investigation of the effects of γ-tACS on working memory in schizophrenia. J Neural Transm.123, 1205–1212. 10.1007/s00702-016-1554-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  HuttA.AtayF. M. (2007). Spontaneous and evoked activity in extended neural populations with gamma-distributed spatial interactions and transmission delay. Chaos Sol. Fract. 32, 547–560. 10.1016/j.chaos.2005.10.091

  • CrossRef
  • Google Scholar

• 30

  HuttA.BuhryL. (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. 10.1007/s10827-014-0512-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  HuttA.BestehornM.WennekersT. (2003). Pattern formation in intracortical neuronal fields. Netw. Comput. Neural Syst. 14, 351–368. 10.1088/0954-898X_14_2_310

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  HuttA.MierauA.LefebvreJ. (2016). Dynamic control of synchronous dynamics in networks of spiking neurons. PLoS ONE11:e0161488. 10.1371/journal.pone.0161488

  • CrossRef
  • Google Scholar

• 33

  JansenB. H.RitV. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybern.73, 357–366. 10.1007/BF00199471

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  KangK.ShelleyM.SompolinskyH. (2003). Mexican hats and pinwheels in visual cortex. Proc. Natl. Acad. Sci. U.S.A. 100, 2848–2853. 10.1073/pnas.0138051100

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  KarK.DuijnhouwerJ.KrekelbergB. (2017). Transcranial alternating current stimulation attenuates neuronal adaptation. J. Neurosci.30, 2266–2216. 10.1523/JNEUROSCI.2266-16.2016

  • CrossRef
  • Google Scholar

• 36

  KlimeschW.SausengP.HanslmayrS. (2007). EEG alpha oscillations: the inhibition-timing hypothesis. Brain Res. Rev.53, 63–88. 10.1016/j.brainresrev.2006.06.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  KunzeT.HunoldA.HaueisenJ.JirsaV.SpieglerA. (2016). Transcranial direct current stimulation changes resting state functional connectivity: a large-scale brain network modeling study. Neuroimage140, 174–187. 10.1016/j.neuroimage.2016.02.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  LafonB.RahmanA.BiksonM.ParraL. C. (2017). Direct Current stimulation alters neuronal input/output function. Brain Stimul.10, 36–45. 10.1016/j.brs.2016.08.014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  LakatosP.KarmosG.MehtaA. D.UlbertI.SchroederC. E. (2008). Entrainment of neuronal oscillations as a mechanism of attentional selection. Science320, 110–113. 10.1126/science.1154735

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  LefebvreJ.HuttA. (2013). Additive noise quenches delay-induced oscillations. Europhys. Lett. 102:60003. 10.1209/0295-5075/102/60003

  • CrossRef
  • Google Scholar

• 41

  LefebvreJ.HuttA.KnebelJ. F.WhittingstallK.MurrayM. M. (2015). Stimulus statistics shape oscillations in nonlinear recurrent neural networks. J. Neurosci. 35, 2895–2903. 10.1523/JNEUROSCI.3609-14.2015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  LindénH.PettersenK. H.EinevollG. T. (2010). Intrinsic dendritic filtering gives low-pass power spectra of local field potentials. J. Comput. Neurosci.29, 423–444. 10.1007/s10827-010-0245-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

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

  • CrossRef
  • Google Scholar

• 44

  Lopes da SilvaF. H.HoeksA.SmitsH.ZetterbergL. H. (1974). Model of brain rhythmic activity. Kybernetik15, 27–37. 10.1007/BF00270757

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  MassiminiM.FerrarelliF.HuberR.EsserS. K.SinghH.TononiG. (2005). Breakdown of cortical effective connectivity during sleep. Science309, 2228–2232. 10.1126/science.1117256

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  MierauA.KlimeschW.LefebvreJ. (2017). State-dependent alpha peak frequency shifts: experimental evidence, potential mechanisms and functional implications. Neuroscience360, 146–154. 10.1016/j.neuroscience.2017.07.037

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  NeulingT.RachS.HerrmannC. S. (2013). Orchestrating neuronal networks: sustained after-effects of transcranial alternating current stimulation depend upon brain states. Front. Human Neurosci. 7:161. 10.3389/fnhum.2013.00161

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  PressW. H.TeukolskyS. A.VetterlingW. T.FlanneryB. P. (2007). Numerical Recipes: The Art of Scientific Computing, 3rd Edn. Cambridge: Cambridge University Press.

  • Google Scholar

• 49

  RobinsonP. A.RennieC. J.WrightJ. J.BahramaliH.GordonE.RoweD. L. (2001). Prediction of electroencephalographic spectra from neurophysiology. Phys. Rev. E63:021903. 10.1103/PhysRevE.63.021903

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  RomeiV.BauerM.BrooksJ. L.EconomidesM.PennyW.ThutG.et al. (2016). Causal evidence that intrinsic beta-frequency is relevant for enhanced signal propagation in the motor system as shown through rhythmic TMS. Neuroimage125, 120–130. 10.1016/j.neuroimage.2015.11.020

  • CrossRef
  • Google Scholar

• 51

  SingerW.GrayC. M. (1995). Visual feature integration and the temporal correlation hypothesis. Annu. Rev. Neurosci.18, 555–586. 10.1146/annurev.ne.18.030195.003011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  SpieglerA.HansenE.BernardC.McIntoshA.JirsaV. (2016). Selective activation of resting-state networks following focal stimulation in a connectome-based network model of the human brain. eNeuro3:ENEURO.0068–16.2016. 10.1523/ENEURO.0068-16.2016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  SpieglerA.KiebelS. J.AtayF. M.KnöscheT. R. (2010). Bifurcation analysis of neural mass models: impact of extrinsic inputs and dendritic time constants. Neuroimage52, 1041–1058. 10.1016/j.neuroimage.2009.12.081

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  SutherlandC.DoironB.LongtinA. (2009). Feedback-induced gain control in stochastic spiking networks. Biol. Cybern.100, 475–489. 10.1007/s00422-009-0298-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  ThutG.MiniussiC. (2009). New insights into rhythmic brain activity from TMS–EEG studies. Trends Cogn. Sci.13, 182–189. 10.1016/j.tics.2009.01.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  UhlhaasP. J.SingerW. (2006). Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron52, 155–168. 10.1016/j.neuron.2006.09.020

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  VarelaF.LachauxJ.-P.RodriguezE.MartinerieJ. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci.2, 229–239. 10.1038/35067550

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  WilsonH. R.CowanJ. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J.12, 1–24. 10.1016/S0006-3495(72)86068-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  WomelsdorfT.ValianteT. A.SahinN. T.MillerK. J.TiesingaP. (2014). Dynamic circuit motifs underlying rhythmic gain control, gating and integration. Nat. Neurosci. 17, 1031–1039. 10.1038/nn.3764

  • Pubmed Abstract
  • CrossRef
  • Google Scholar