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

$\begin{array}{ll}{X}_n^j\left(t\right)\to Poisson\left(f\left[u_n^j\left(t\right)\right]\right)& (1)\end{array}$

where $X_n^j(t)={\displaystyle \sum _{t_l}{\delta }_n^j}(t-t_l)$ is the ensemble spike train of the jth patch and n = e, i indicate excitatory and inhibitory populations, respectively. The firing rate 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\left[u_n^j\right]={\left(1+exp\left[-\beta \left(u_n^j-h\right)\right]\right)}^{-1}$ i.e., the firing rate probability approaches f = 1 for large membrane potentials, where the gain is β > 0 and the firing rate threshold is h. This defined, the model combines both the spiking of single cells as well as a dependence on the firing rate of the whole population. Such a hybrid cortical model thus combines both spiking and rate driven dynamics. The excitatory and inhibitory potentials $u_e^j\left(t\right)$ and $u_i^j\left(t\right)$ represent ensemble-averaged potentials proportional to averaged dendritic currents. They obey the set of non-linear stochastic equations

$\begin{array}{ll}{\alpha }_n^{-1}\frac{d{u}_n^j\left(t\right)}{dt}=L\left[u_n^j\left(t\right)\right]+{\displaystyle \sum _m}{G}_{nm}^j\left(t\right)+\sqrt{2D_n}{\xi }_n^j\left(t\right)+S\left(t\right)& (2)\end{array}$

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

$\begin{array}{ll}{G}_{nm}^j\left(t\right)={\displaystyle \sum _{k=1}^{N_m}}{W}_{nm}^{jk}\left(c\right)·PS{P}_m^k\left(t-{\tau }^{jk}\right)& (3)\end{array}$

where $PS{P}_m^k\left(t\right)$ refers to mean post-synaptic potential of patch k in population m at time t. Interactions between subnetworks k and l are subjected to intracortical propagation delaysτ^jk = |x(k)−x(i)| v⁻¹, 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

$\begin{array}{ll}PS{P}_m^k\left(t\right)={\displaystyle \underset{o}{\stackrel{t}{\int }}}{X}_m^k\left(s\right)\frac{1}{a_m}{e}^{-\frac{t-s}{a_m}}ds& (4)\end{array}$

with synaptic time constant a_m.

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

$\begin{array}{ll}{W}_{nm}^{jk}\left(c\right)=w_{nm}^o\left(c\right)exp\left[-{\sigma }_{n,m}^2|x\left(j\right)-x\left(k\right)|\right]& (5)\end{array}$

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

$\begin{array}{ll}A\left(t\right)=\frac{1}{N_e}{\displaystyle \sum _{k=1}^{N_e}}{\varphi }_e^k{u}_e^k\left(t\right)+\frac{1}{N_i}{\displaystyle \sum _{k=1}^{N_i}}{\varphi }_i^k{u}_i^k\left(t\right)& (6)\end{array}$

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

TABLE 1

Table 1. Cortical microcircuit network model parameters.

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

$\begin{array}{ll}{G}_{ee}^j,G_{ei}^j\ll G_{ie}^j,G_{ii}^j.& (7)\end{array}$

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

$\begin{array}{l}{\alpha }_e^{-1}\frac{d{u}_e^j\left(t\right)}{dt}\approx L\left[u_e^j\left(t\right)\right]+G_{ie}^j\left(t\right)+\sqrt{2D_e}{\xi }_e^j\left(t\right)+S\left(t\right)\\ {\alpha }_i^{-1}\frac{d{u}_i^j\left(t\right)}{dt}\approx L\left[u_i^j\left(t\right)\right]+G_{ii}^j\left(t\right)+\sqrt{2D_i}{\xi }_i^j\left(t\right)+S\left(t\right)& (8)\end{array}$

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

$\begin{array}{l}{G}_{ii}^j\left(t\right)\approx {\sum }_{k=1}^{N_m}{W}_{ii}^{jk}\left(c\right)·f\left[u_i^k\left(t-{\tau }^{jk}\right)\right]\\ \approx w_{ii}^o{\sum }_{k=1}^{N_m}f\left[u_i^k\left(t-\bar{\tau }\right)\right]& (9)\end{array}$

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

$\begin{array}{ll}{\alpha }^{-1}\frac{d}{dt}{U}^j\left(t\right)=L\left[U^j\left(t\right)\right]+g{N}^{-1}{\sum }_{j=1}^{N}f\left[U^j\left(t-\bar{\tau }\right)\right]+S\left(t\right)& (10)\end{array}$

where ${\alpha }^{-1}\equiv {\alpha }_i^{-1}$ and $U^j\left(t\right)\equiv u_i^j\left(t\right)$ was introduced to distinguish the more detailed cortical microcircuit model and the reduced model in the subsequent calculations. The mean synaptic connectivity $w_{ii}^o=g{N}^{-1}$ has also been introduced to indicate the average evaluated over all possible pairs of inhibitory neurons. According to the derivations above the mean synaptic action is inhibitory with g < 0. The non-linear response function above remains the same with f[u] = (1+exp[−β(u−h)])⁻¹. 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.,

$\begin{array}{ll}{U}^j\left(t\right)=Ū\left(t\right)+V^j\left(t\right)& (11)\end{array}$

where Ū is given by

$\begin{array}{ll}Ū\left(t\right)=N^{-1}{\displaystyle \sum _{i=1}^N}{U}^j\left(t\right)\equiv <U{>}_N& (12)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}{V}^j=-V^j+S\left(t\right)-{\mu }_S,& (13)\end{array}$

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

$\begin{array}{ll}{\mu }_S={\displaystyle \underset{t}{\stackrel{t+T}{\int }}}S\left(s\right)ds& (14)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}Ū\left(t\right)=-Ū\left(t\right)+gF\left[Ū\left(t-\bar{\tau }\right)\right]+{\mu }_S.& (15)\end{array}$

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

$\begin{array}{ll}F\left[Ū\right]={\displaystyle \underset{\Omega \left(\nu \right)}{\stackrel{}{\int }}}f\left[Ū+V\right]\rho \left(V\right)dV& (16)\end{array}$

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

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

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

FIGURE 2

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

FIGURE 3

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

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

Impact of Stimulus Waveform on Alpha Oscillations: Theoretical Insights

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

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

FIGURE 4

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

TABLE 2

Table 2. Reduced model parameters.

1. Pulse Train Stimulation

Let us first consider the stimulus waveform

$\begin{array}{ll}S\left(t\right)=S{\sum }_n\delta \left(t-\frac{2\pi }{r}n\right),& (17)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}{V}^j=-V^j+S{\sum }_n\delta \left(t-\frac{2\pi }{r}n\right)-{\mu }_S& (18)\end{array}$

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)

$\rho \left(V\right)=\{\begin{array}{c} \frac{r}{V-{\mu }_S} | V\in [S{e}^{-\frac{1}{r}},S]\\ 0 | \text{otherwise}& (19)\end{array}$

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

$\begin{array}{l}F\left[Ū\right]=r{\displaystyle \underset{S}{\stackrel{S{e}^{-\frac{1}{r}}}{\int }}}\frac{H\left[Ū+V-h\right]}{V}dV\\ =-rH\left[Ū+S-h\right]\left(ln\left[-Ū\right]-ln\left[S\right]\right)+rH\left[Ū+Sexp\left(-\frac{1}{r}\right)-h\right]\\ \left(ln\left[-Ū\right]-ln\left[Sexp\left(-\frac{1}{r}\right)\right]\right).& (20)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}Ū\left(t\right)=-Ū\left(t\right)+gF\left[Ū\left(t-\bar{\tau }\right)\right]+Sr.& (21)\end{array}$

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

$\begin{array}{ll}S\left(t\right)=Ssin\left(2\pi {\omega }_{s}t\right)& (22)\end{array}$

with angle frequency ω_s. Fluctuations around the mean obey

$\begin{array}{ll}\frac{d}{dt}{V}^j=-V^j+Ssin\left(2\pi {\omega }_{s}t\right)& (23)\end{array}$

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)

$\begin{array}{ll}\rho \left(V\right)=\frac{1}{\pi \sqrt{{\left(S/2\pi {\omega }_s\right)}^2-V^2}}& (24)\end{array}$

and consequently the effective non-linearity is given by

$\begin{array}{l}\tilde{F}\left[Ū\right]\approx {\displaystyle \underset{-\mu }{\stackrel{\mu }{\int }}}\frac{H\left[Ū+V-h\right]}{\pi \sqrt{{\left(\frac{S}{2\pi {\omega }_s}\right)}^2-V^2}}dv=\frac{1}{\pi }{sin}^{-1}\left(2\pi {\omega }_SŪ/S\right)\\ +\frac{1}{2},-1\le 2\pi {\omega }_SŪ/S\le 1\\ \tilde{F}\left[Ū\right]=1,2\pi {\omega }_SŪ/S>1\\ \tilde{F}\left[Ū\right]=0,2\pi {\omega }_SŪ/S<-1& (25)\end{array}$

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

$\begin{array}{ll}S\left(t\right)=S_j\left(t\right)=\sqrt{2D}{\xi }_j\left(t\right),& (26)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}{V}^j=-V^j+\sqrt{2D}{\xi }_j\left(t\right)& (27)\end{array}$

where μ_S = 0. The associated probability density function ρ(V) is a symmetric Gaussian

$\begin{array}{ll}\rho \left(V\right)=\frac{1}{\sqrt{2\pi D}}exp\left[-\frac{V^2}{2D}\right].& (28)\end{array}$

Convolving with the threshold non-linearity yields

$\begin{array}{ll}\tilde{F}\left[Ū\right]\approx \frac{1}{\sqrt{2\pi D}}{\displaystyle \underset{-\mu }{\stackrel{\mu }{\int }}}H\left[Ū+V-h\right]exp\left[-\frac{V^2}{2D}\right]dv=\frac{1}{2}+\frac{1}{2}\text{erf}\left[\frac{Ū}{\sqrt{2D}}\right].& (29)\end{array}$

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

$\begin{array}{ll}S\left(t\right)=S=constant& (30)\end{array}$

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

$\begin{array}{ll}F\left[U\right]=H\left[U-h\right]\approx f\left[U\right]& (31)\end{array}$

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

$\begin{array}{ll}Ū\to Ū-S& (32)\end{array}$

which leads to the mean–field dynamics

$\begin{array}{ll}\frac{d}{dt}Ū\left(t\right)=-Ū\left(t\right)+g\tilde{F}\left[Ū\left(t-\bar{\tau }\right)\right]& (33)\end{array}$

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

Impact on Alpha Peak Frequency

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

$\begin{array}{ll}{Ū}_o=gF\left[{Ū}_o\right]+{\mu }_S.& (34)\end{array}$

The linearized dynamics around the steady state ū_o is

$\begin{array}{ll}\frac{d}{dt}Ū\left(t\right)=-Ū\left(t\right)+R\left[{Ū}_o\right]Ū\left(t-\bar{\tau }\right)& (35)\end{array}$

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

$\begin{array}{ll}J={\displaystyle \underset{0}{\overset{2\pi /\omega }{\int }}}\left(\frac{d}{dt}Ū-F\left[Ū\left(t-\bar{\tau }\right)\right]\right)·\mathrm{\text{cos}}\left(\omega t\right)dt& (36)\end{array}$

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

$\begin{array}{ll}{J}_{lin}={\displaystyle \underset{0}{\overset{2\pi /\omega }{\int }}}\left(\frac{d}{dt}Ū-R\left[{Ū}_o\right]Ū\left(t-\bar{\tau }\right)\right)·\mathrm{\text{cos}}\left(\omega t\right)dt& (37)\end{array}$

Using the ansatz Ū(t) = Acos(ωt)+Ū_o, one obtains for the first iteration

$\begin{array}{ll}{J}_{lin}=-\frac{A\pi \left(R\left[{Ū}_o\right]cos\left(\omega \tau \right)-1\right)}{\omega }& (38)\end{array}$

Setting J_lin = 0 and solving for ω yields an approximation of the linear frequency as a function of R[ū_o] and delay τ i.e.

$\begin{array}{ll}\omega \approx \frac{arccos\left(\frac{1}{R\left[{Ū}_o\right]}\right)}{\bar{\tau }}& (39)\end{array}$

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., $\bar{\tau }=25$ms). Although R[Ū_o] cannot always be computed analytically due to the implicit condition for the equilibrium (34), some specific cases such as Gaussian white noise for instance, remain tractable and accurate (Hutt et al., 2016), and can otherwise be computed numerically.

FIGURE 5

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

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

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

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

FIGURE 6

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

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

FIGURE 7

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

Discussion

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

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

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

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

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

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

Author Contributions

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

Conflict of Interest Statement

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

Acknowledgments

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

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Womelsdorf, T., Valiante, T. A., Sahin, N. T., Miller, K. J., and Tiesinga, P. (2014). Dynamic circuit motifs underlying rhythmic gain control, gating and integration. Nat. Neurosci. 17, 1031–1039. doi: 10.1038/nn.3764

PubMed Abstract | CrossRef Full Text | Google Scholar