Dynamic Control of Synchronous Activity in Networks of Spiking Neurons

Oscillatory brain activity is believed to play a central role in neural coding. Accumulating evidence shows that features of these oscillations are highly dynamic: power, frequency and phase fluctuate alongside changes in behavior and task demands. The role and mechanism supporting this variability is however poorly understood. We here analyze a network of recurrently connected spiking neurons with time delay displaying stable synchronous dynamics. Using mean-field and stability analyses, we investigate the influence of dynamic inputs on the frequency of firing rate oscillations. We show that afferent noise, mimicking inputs to the neurons, causes smoothing of the system’s response function, displacing equilibria and altering the stability of oscillatory states. Our analysis further shows that these noise-induced changes cause a shift of the peak frequency of synchronous oscillations that scales with input intensity, leading the network towards critical states. We lastly discuss the extension of these principles to periodic stimulation, in which externally applied driving signals can trigger analogous phenomena. Our results reveal one possible mechanism involved in shaping oscillatory activity in the brain and associated control principles.


Introduction
Brain signals are rife with oscillatory spectral patterns. These rhythmic features, uncovered through both intracranial and non-invasive recordings, have been shown to correlate strongly with cognitive processes, memory, and sensorimotor behavior [1,2,3], and are thus believed to be dynamic signatures of specific neural computations. As such, oscillatory activity is ubiquitous throughout the nervous system and represents the focus of an effervescent area of research [4,5,6].
Brain oscillations are however far from being static. Indeed, cortical rhythms are commonly subjected to sudden shifts induced by variations in behavior and cognitive states. Such spectral transitions are notably observed across normal sleep stages [7] or during the recruitment of attention [8]. The variability of oscillatory neural activity within the gamma band has been thoroughly studied and linked to changes in visual stimuli statistics and timely adjustments in local synaptic wiring [9,10,11]. However, shifts in brain oscillatory activity are also reliably observed at slower frequencies. Alpha oscillatory activity, in particular, has been found to be highly volatile [12]. It has been found that the iAPF accelerates during cognitive [13], memory [14] and sensorimotor [15] task performance as well as following a strenuous bout of physical exercise [16]. In addition, recent studies provide strong evidence that alpha oscillations are one candidate mechanism for gating the temporal window of sensory integration, and thus dictating the resolution of conscious sensory updating. Specifically, deliberate alterations of the iAPF within individual subjects, induced by transcranial alternating current stimulation, has been found to influence visual stimuli [17]. Consistent with this, individuals with higher iAPF have vision with finer temporal resolution, and within an individual, spontaneous fluctuations in iAPF predict visual perception [18]. Furthermore, in a temporal cueing task, forming predictions about when a stimulus will appear can instantaneously bias the phase of ongoing alphaband oscillations toward an optimal phase for stimulus discrimination [19]. Taken together, these findings suggest that the magnitude of the iAPF is indicative of the level of arousal/attention, preparedness and performance of recruited cortical nets, in which "the faster the better". Given that alpha activity operates on much broader spatial and temporal scales [20], spectral transitions observed within those frequency ranges likely relies on more global and distributed mechanisms.
To this day, the mechanism supporting these larger scale frequency transitions has been poorly understood. A key question is whether such transitions could be triggered by external stimulation. Recent studies have indeed shown that weak electric fields can perturb individual alpha oscillations and have a direct effect on visual stimulus perception [17,21] and task performance by reinforcing endogenous slow-wave rhythms [1,22].
To explore this question, we here investigate a non-linear network of spiking neurons with time delay, exhibiting alpha-like oscillatory activity. Our results reveal that despite the noisiness of the connectivity, stimuli implement an online gain-control mechanism where the peak frequency reflects the activation state of the neurons. We show that neural inputs, here modeled by noise, change the shape of the neuron response function, significantly changing the systems equilibria and stability. Using mean-field analysis of the network collective dynamics, we demonstrate that noise causes the system to shift from slow non-linear oscillations to fast linear oscillations, bringing the system towards a critical state. We also derive a frequency tuning curve that relates the network's synchronous frequency to the noise intensity driving its constituent neurons. We lastly explore how these results also apply to periodic stimuli, opening new perspectives on how external brain stimulation can be used to control neural synchronous activity.

Model
In the present work, we analyze the dynamics of a generic network of spiking neurons (Fig 1A) whose membrane potential u i (t) evolves according to the following set of non-linear differential equations where α is the membrane time constant, w ij = [W] ij are synaptic weights and where τ is a mean conduction delay. The membrane potential u i (t) represents the deviation from the neurons resting potential that is present in the absence of synaptic and external input. The presynaptic spike trains X j ðtÞ ¼ S ft l g dðt À t l Þ obey the non-homogenous Poisson processes X i ! Poisson (f[u i ]) with rate f and the Dirac distribution δ(t). The firing rate function f has a non-linear sigmoid shape and is defined by f[u i ] = (1 + exp[−βu i ]) −1 ,i.e. the maximum firing rate approaches f = 1 for large membrane potentials. All neurons are subjected to afferent pre-synaptic inputs. The synaptic connectivity scheme was randomly set (Fig 1B), such that where g is the mean synaptic strength, s is the weight variance and n ij are zero-mean independent Gaussian white noise such that < n ij n kl > NxN = δ ik δ jl with the Kronecker symbol δ nm and where < > NxN is an average evaluated over all possible pairs of indices of the matrix W. We also consider other sources of synaptic inputs which we model as stochastic elements x i ðtÞ given the independent Gaussian white noise processes ξ i with zero mean and < ξ i ξ j > T = δ ij , where < > T is an average evaluated over an epoch of duration T or over an ensemble of realizations of processes. Such noise is meant to represent the effect of synaptic bombardment on neural membrane potentials. Such inputs can be shown to be well approximated by Gaussian processes in the diffusion limit case [23,24] and this is the approach we use in the following analysis. In the present work, we also consider a network of neurons whose spatial mean synaptic action is inhibitory with < 0.
Under the choice of a specific synaptic connectivity and other model parameters, the network spontaneously engages in synchronous activity at a baseline frequency of 10Hz while external input is absent. However, this frequency is found to be highly volatile in the presence of stochastic input: if the noise level is increased, the network spiking activity becomes more irregular and the synchronous frequency increases, cf. Fig 1C. The power spectrum of the network mean activity is plotted in Fig 1D, where peaks can be observed at the systems natural frequency and higher harmonics. As noise intensity is increased in the system, the peak frequency and associated harmonics gradually shift towards higher a frequency range.

Mean Field Representation for Stochastic Input
To better understand the impact of inputs on the dynamics, let us derive mean-field equations for our system. In the limit where the number of neurons is large, i.e. N ! 1, and the mean firing time scale 1/f is much smaller than the time scale of dendritic currents [25] we can use the response function f to approximate the rate at which recurrent pre-synaptic inputs perturb the activity of a neuron. To this end, one averages u i (t) over a very short time window and introduces a so-called coarsening in time [25,26]. This well-established transformation allows to translate population spiking activity to population rate dynamics In the following, we use the same symbol for the original and the temporally coarse-grained membrane potential for notational simplicity. Let us further assume that emerging oscillations occur in a mean-driven regime in which the local dynamics can be seen as small independent fluctuations around the network mean activity i.e.
where the network mean activity is given by and < > N is an average performed over the N unit of the network. In the following, we re-scale time by αt ! t for notational simplicity. As an ansatz, local fluctuations v i from the mean obey the Ornstein-Uhlenbeck processes Then, using Eq (4) above, and taking the mean over N neurons Then, as N ! 1 [27], where ρ is the probability density function of the solution of Eq (6), i.e. a zero-mean Gaussian distribution with variance var for the expectation value E of a product of two statistically independent random variables X and Y. This applies in Eq (8) Combining previous results, the mean-field dynamics of our spiking network can be well approximated by the scalar non-linear delay-differential equation Eq (10) shows that pre-synaptic noise generically results in a linearization of the neurons response function [28,29] cf. Fig 2, bottom panels. This linearization shapes the input-output relationship of driven neurons, but also alters the features displayed by emergent activity patterns such as synchronous oscillations [30].

Stochastic Stability Analysis
To investigate the frequency tuning observed in Fig 2, we take a closer look at the solutions of Eq (10). Synchronous activity in our network emerges as a consequence of an expected supercritical Hopf bifurcation commonly seen in recurrent delayed nets [31]. The single fixed point " where the last equation assumes small noise intensity D. Fig 3A shows that the mean-field fixed point " u o decreases with increasing noise level. We use the knowledge about the noise-dependent fixed point " u o to better understand the network stability. The linearization of the mean-field Eq (10) about the fixed point " u o from Eq (11) yields with deviations w from the fixed point where R½" u o ; D ¼ " bility. Stability of the network equilibrium " u o is determined by setting wðtÞ ¼ũe lt with l 2 C leading to the characteristic equation Solutions to Eq (13) define the spectrum of the linearized system in Eq (12). Collective oscillatory solutions form in the network if the susceptibility reaches a critical value R c for λ = i ω c satisfying the complex relationship where o c ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi R c À 1 p is the critical Hopf frequency i.e. the frequency of network synchronous oscillations close to the instability.
As seen in Fig 3B, in the absence of noise, the susceptibility is well below the critical value, and the network exhibits strong non-linear oscillations. This implies that R < R c for D = 0 and the system evolves in a nonlinear limit cycle oscillation beyond a supercritical Hopf bifurcation while the fixed point is asymptotically unstable. As such, in the weak noise limit, the system is set robustly in the synchronous state.  However, upon the presence of noise, linearization of the neurons response function translates into a gradual increase of the susceptibility towards R c : global synchronous oscillations in the network not only accelerate under the action of noise, but also becomes more linear. As such, noise brings the system closer to the bifurcation threshold and hence moves the network towards a critical asynchronous state. The effect of noise on the characteristic eigenvalue spectrum is plotted in Fig 3C. As noise intensity increases, the eigenvalues undergo a gradual shift towards the left hand side of the complex plane, bringing pairs of eigenvalues closer to the imaginary axis and thus closer to the bifurcation threshold. This occurs because the system's susceptibility approaches the critical susceptibility under the effect of noise. As such, afferent inputs engage the network and drive it towards the asynchronous state and the transition trajectory in parameter space is characterized by an gradual increase in the network peak frequency.

Stochastic Frequency Tuning
Finding explicit frequency relationships in the fully non-linear regime is challenging. Yet, to approximate the dependence of the network frequency on the noise driving the neurons, one might use the Galerkin method [32,33]. Whenever f[u] % H[u] holds for sufficiently large values of the gain β, this approach seeks to find a frequency ω minimizing the measure where F½" u ¼ À " uðtÞ þ "  For J = 0, one can solve Eq (16) for ω > 0 to obtain a first order approximation for the stochastic frequency Expanding to third order for weak noise, i.e. D % 0 leading to the frequency tuning curve where Δ(D) is a noise-induced shift with Δ(D) > 0. Fig 4 shows the frequency with respect to the noise intensity. Together with the results above, while approximate, demonstrate that under the action of noise, the network oscillation frequency gradually shifts from a non-linear baseline frequency (18) and (19) thus imply that as noise increases in the system, the network transits from slow non-linear rhythms to fast linear oscillations. This can also be understood by looking at the susceptibility, which gauges the relative influence of recurrent interactions in the network, plotted in Fig 3B. As susceptibility decreases under the action of noise, the system accelerates and shifts from a recurrent deeply synchronous state to an asynchronous input driven regime.

Mean Field Analysis for Periodic Stimulation
We have seen that noise shapes the stability and oscillatory features of non-linear networks by linearizing the neurons response function. However, this feature is not exclusive to noisy inputs. Indeed, numerous studies have shown that, in addition to changes in brain state, ongoing oscillatory activity can be modulated by noninvasive stimulation, something that is increasingly capitalized upon in basic research and clinical practice [34,35,36,37,38,39]. Recent results have further shown that exogenous electric rhythmic stimulation (i.e. periodic forcing), in addition to resonance and entrainment, can also provoke non-linear acceleration of endogenous oscillations and shift the baseline frequency of driven neural systems [40], causing resonance curves and Arnold tongues to bend in stimulation parameter space ( [40] cf Fig 4). To explore the mechanism behind this non-linear effect, we here revisit the analysis performed in the stochastic case and adapt it to the presence of periodic forcing.
In the presence of global periodic stimulation with frequency f s and amplitude I 0 , the network dynamics obeys First, we consider activity course-grained in time with similar to the stochastic case. Then we assume two temporal scales in the evolution of the potentials with u j = m j + v j : the slow mode m j and the fast mode v j . Inserting this relation into Eq (20), we obtain and d dt m i ðtÞ ¼ À m i ðtÞ þ N À 1 S N j¼1 w ij f ½m j ðt À tÞ þ v j ðt À tÞ ð23Þ Since the experimental observation reflects the average activity in the neural ensemble, we consider spatially homogeneous slow activity m j ¼ " uðtÞ. Then averaging over the time interval of one short stimulus cycle, one can express the spatio-temporal mean dynamics " uðtÞ in the adiabatic regime for fast stimuli Here <Á> T,N denotes a spatial average and a time average taken over a time interval T small compared to the dynamics of the network i.e. for 1/f s T ( 2π/ω o , i.e. for the local activity u i (t). Again, this implies that the driving frequency f s is large compared to the systems self-sustained oscillation frequency. By virtue of this time scale separation, it is reasonable to assume stationarity in the time interval of duration T. The solutions of Eq (22) for large frequencies and large times t ! 1 obey v i ðtÞ % À I 0 2pf s cosð2pf s tÞ; ð26Þ i.e. the fluctuations about the spatial mean synchronize and converge to a single solution. Then with the probability density with μ = I 0 /2πf s . Here, we have used the time scale separation between the slow evolution of " uðtÞ and the fast synchronous fluctuations. Eq (28) has the same form as Eq (8) in the case of stochastic stimulation. Again, similar to the stochastic case, the external stimulation leads to a convolution of the nonlinear response function f with the probability density function of the stimulus-induced fluctuations about the spatial mean. We note that, in contrast to the stochastic case, the notion of a probability density function in the present deterministic system may appear counter-intuitive. This interpretation however is reasonable since ρ(v) is proportional to the residence time of the oscillation v i (t) at amplitude v [41] that motivates the interpretation of a probability density.

Stability and Frequency Tuning
Taking a closer look at the slope of the new response functionF 0 ½u ¼ dF=du, we find for À m " u m, and dF0 dm ¼ 0 otherwise. Hence the nonlinear response function flattens and becomes increasingly linear for increasing μ similar to the stochastic case for increased noise level. Fig 5 (bottom panels) showsF for three different stimulus amplitudes, i.e. different values of m, confirming this analytical finding.
To gain insight into the dynamics of " u, we consider the fixed point " u o and small deviations about it. Similar to Eq (12), the corresponding characteristic roots are defined by Since dR dm > 0, similar to the stochastic case fast external periodic driving moves the system from a nonlinear regime far from the bifurcation threshold to a linear regime close to to the stability threshold, cf. Fig 3C. Using the Galerkin approach as in Eq (15), one obtains an analogous expression to Eq (17) for the oscillation frequency but for the periodic forcing case for small periodic driving amplitudes I 0 and with the fixed point " u o ¼F½" u o % pm " w=2ðmp À " wÞ for small μ. Then after few calculus steps one finds for small values of μ. Hence linearising the response function by fast external stimulation increases the oscillation frequency of the system. Fig 5 shows numerical simulations of Eq (30) for three different driving amplitudes. The oscillation frequency increases with increasing driving amplitude in accordance to the analytical result in Eq (28). This finding resembles the stochastic case for increasing noise level. To summarize the results above, periodic forcing not only interacts with the dynamics of the network by resonance or entrainment, but also shapes its activity via non-linear effects. This thus implies that the use of weak, high frequency stimulation can mediate frequency transitions in a similar fashion as stochastic inputs.

Discussion
Previous experimental and theoretical work has shown that gamma-like oscillatory activity is highly dynamic, changing according to stimulus intensity [42]11), spatial features [43,44] and is strongly correlated to the phase of slower frequencies [45,46]. Such gamma oscillations have been shown to build on highly local circuits shaping the timing of interactions between excitatory pyramidal cells and inhibitory interneurons [11]. Slower frequencies, such as alpha activity, have also been shown to be variable. Shifts in the peak alpha frequency have been reliably reported during changes in attentional states [13,18,19], during sensorimotor task performance [15], and following intense physical exercise [16].
Alpha oscillations have been shown to engage more spatially extended connections [47], in which propagation delays play an important role [20], suggesting that the mechanisms involved in shaping the peak alpha frequency is different from the one involved for faster, more local frequencies such as gamma for example. Most of the research on synchronous neural dynamics has been devoted to the study of input-induced transitions in-and out of synchronous states, where network firing rate oscillations emerge in presence of strongly correlated drive [30,48,49]. In contrast, we here characterize a smooth approach towards a delayed-induced bifurcation in which the equilibria, stability and peak frequency are impacted by noise or external periodic driving. Using mean-field approaches, we have detailed the stability of the network oscillatory states and derived a frequency tuning curve that relates the frequency of synchronous oscillations to the intensity of the input driving the neurons.
While we have considered constant values of the input intensity in our analysis, results extend to time-varying inputs as well. For instance, the level of noise D(t) driving neurons would vary according to particular sets of stimuli and\or tasks. In such a case, fluctuations in the noise intensity would be mirrored by concomitant changes in the network peak frequency. What are implications of this dynamic behavior with respect to neural coding? Similarly, experimental setups involving electric periodic stimulation, such as in deep brain stimulation [50,51,52], electric stimulation [53,54].or visual stimulation [55], will induce a shift of frequency [40].
According to the framework we detailed, shifts in the magnitude of inputs to the neurons translate into changes in synchronous oscillations frequencies. Network mean output activity can thus be described by the heuristic relationship where f o is the baseline frequency, Δ is an input-dependent frequency shift, A is an amplitude parameter (which may depend on input parameters) and ϕ is a random phase. From this perspective, the frequency tuning mechanism can be seen as implementing a frequency modulation (FM) communication scheme coding for input intensity. This hypothesis has been considered on numerous occasions [56,30], and lately gained momentum alongside the development of the "Communication Through Coherence" (CTC) hypothesis in which neural assemblies engage into long-distance communication based on the coherence between oscillatory states [57]. Hence, we hypothesize that the stochastic frequency tuning explored here could potentially implement a gating mechanism allowing the routing of information towards different distal networks, and thus represent an aspect of large-scale neural coding.