## Abstract

Rhythmic activity plays a central role in neural computations and brain functions ranging from homeostasis to attention, as well as in neurological and neuropsychiatric disorders. Despite this pervasiveness, little is known about the mechanisms whereby the frequency and power of oscillatory activity are modulated, and how they reflect the inputs received by neurons. Numerous studies have reported input-dependent fluctuations in peak frequency and power (as well as couplings across these features). However, it remains unresolved what mediates these spectral shifts among neural populations. Extending previous findings regarding stochastic nonlinear systems and experimental observations, we provide analytical insights regarding oscillatory responses of neural populations to stimulation from either endogenous or exogenous origins. Using a deceptively simple yet sparse and randomly connected network of neurons, we show how spiking inputs can reliably modulate the peak frequency and power expressed by synchronous neural populations without any changes in circuitry. Our results reveal that a generic, non-nonlinear and input-induced mechanism can robustly mediate these spectral fluctuations, and thus provide a framework in which inputs to the neurons bidirectionally regulate both the frequency and power expressed by synchronous populations. Theoretical and computational analysis of the ensuing spectral fluctuations was found to reflect the underlying dynamics of the input stimuli driving the neurons. Our results provide insights regarding a generic mechanism supporting spectral transitions observed across cortical networks and spanning multiple frequency bands.

## Introduction

Oscillations are an information-rich mechanism of neural signaling (Engel et al., 2001; Varela et al., 2001; van Atteveldt et al., 2014). Rhythmic neural activity involves intermittent synchronous patterns between distant yet simultaneously activated regions (Wang, 2010) and have been suggested to shape communication timing between neural populations (Fries, 2005) and to be centrally involved in the binding of sensory representations (Engel and Singer, 2001; Engel et al., 2001). However, cortical networks exhibiting synchronous dynamics cannot be considered in isolation; neural populations are continually subjected to myriad (stochastic) influences. Afferent fluctuating stimuli impact the dynamics of neural networks both by recruiting individual neurons, and also by shaping collective synchronous activity of neuronal populations.

The nonstationarity of power spectra in the brain has been robustly observed across neurophysiologic measures both intracranially (Schroeder and Lakatos, 2009; Whittingstall and Logothetis, 2009; Ray and Maunsell, 2010; Musall et al., 2014) and noninvasively at the scalp surface (Van Zaen et al., 2010, 2013; Thut et al., 2012). The nature and features of the input shape the peak frequency and power of ongoing cyclic activity among neurons in a nontrivial and time-dependent way, operating on multiple frequency bands (Cohen, 2014). Fluctuating coherently with stimulation dynamics, the frequency of ongoing neural oscillations typically increases alongside input intensity (Whittington et al., 1995; Ermentrout and Kopell, 1998). Yet, the mechanisms responsible for this phenomenon are currently poorly understood.

Changes in the peak frequency displayed by synchronous neural ensembles are typically ascribed to alterations in circuit features, such as network structure or synaptic timescales. However, the occurrence of such changes over timescales as short as 100 ms suggests that other input-dependent mechanisms may underlie rapid fluctuations in power spectra, rather than neuromodulatory and/or plastic changes in circuitry. Mechanisms involved in the generation and shaping of gamma-type activity have been well characterized (Whittington et al., 1995; Wang and Buzsáki, 1996; Ray and Maunsell, 2010) and are considered a signature of highly local synaptic processing (Jadi and Sejnowski, 2014). However, these mechanisms are likely distinct from those supporting modulations of slower oscillations (Lakatos et al., 2008; Schroeder and Lakatos, 2009; Burgess, 2012). Alpha activity is believed to rely on larger scale processes (Hindriks et al., 2014) supported by delayed network interactions (Cabral et al., 2014) operating over longer timescales (Haegens et al., 2014). Here, we specifically investigate how such delay-induced alpha activity is modulated by changes in stimulation dynamics and demonstrate how peak frequency and power reflect the activation state of the neurons on a larger spatial scale. To do this, we investigated mathematically the spectral response of a simple, yet nonlinear neural population model subjected to stochastic spike-like inputs. Using mean-field analysis, our results reveal a dependence of the system's effective nonlinearity on stimulus statistics. Specifically, they were shaped by changes in inputs while all other parameters were kept constant, such that alpha-type oscillatory activity was fully regulated by the dynamics of the stimulus. Time-dependent fluctuations in peak frequency and power can therefore result from a fully generic property of synchronous nonlinear networks, while echoing the dynamics and statistics of underlying inputs.

## Materials and Methods

#### Model description

To analyze the spectral response of synchronous networks to inputs, we deliberately chose a simple, yet nonlinear, sparse and random system to enforce the generality of the results. Using this model, described in detail below, we developed a mean-field representation of the dynamics to determine the effect of stimulus statistics on the activity pattern of the modeled neurons. We then performed a thorough stochastic stability analysis while the system remained close to oscillatory instabilities. Throughout, we studied the effect of changes in input statistics, and no other parameter was varied.

In what follows, we deliberately consider a simplified network of synchronous neural units subjected to delayed and nonlinear reciprocal interactions and stimulated by another, upstream population via independent trains of modeled action potentials. The associated representation of our model is aimed to provide a description of brain dynamics at the mesoscopic (centimeter) scale (Freeman, 1975; Wright and Liley, 1995; Jirsa and Haken, 1997), in which slower rhythmic activity modulates faster, more local oscillations (Osipova et al., 2008). A schematic illustration is shown in Figure 1. The quantity *u _{i}*(

*t*) with

*i*= 1,…,

*N*represents the mean somatic membrane activity of the

*i*th subunit in the network, and obeys the nonautonomous set of dimensionless evolution equations: where α = 50 ms defines the mean synaptic rate of the neuronal population,

**(**

*I**t*) is a vector of stochastic inputs of dimension

*N*, and the vector

**of dimension**

*f**N*denotes the nonlinear response functions of the neurons. The synaptic efficacies of the neurons are identical and given by

*g*=

*g*

_{o}N^{−1}, where again

*N*is the number of units in the network. Without constraining generality the nonlinear response function displays the sigmoid-like shape

*[*

**f***u*]

*=*

_{j}*f*[1 + exp(−β

_{o}*u*)]

_{j}^{−1}, where β = 100 is the steepness, or gain, of the neural responses. The quantity

*f*, here set to 100 Hz, is the maximum firing rate of the neurons. For simplicity, network interactions are conditioned by a unique and constant propagation delay τ = 30 ms, which approximates well the distribution of interaction latencies in sparse and random networks (Roxin et al., 2005). This approximation is, however, not mandatory, and τ can be thought to represent the maximal (i.e., longest) interaction latency. Units in the network are also subjected to independent time-varying stimuli. Spiking activity of individual units in the network follows a nonhomogeneous Poisson process, i.e.,

_{o}*X*(

_{i}*t*)→

*Poisson*(

*f*(

*u*(

_{i}*t*))).

The network's synaptic connectivity matrix ** W**(

*c*) is depicted in Figure 1

*B*. Its structure is assumed to be sparse and random with connection probability

*c*= 0.8. Random synaptic weights, sampled from a density ρ

*(*

_{W}**), are assumed to be locally excitatory and distally inhibitory: synaptic weights [**

*W***]**

*W**=*

_{ij}*w*for which |

_{ij}*i*−

*j*| < r, given some radius

*r*, are made positive valued (within the interval [0,1]), while all others are negative (within the interval [−1, 0]). Synaptic weights are uniformly distributed over their respective intervals, i.e., ρ

*(*

_{W}*w*) = 1. In addition, an amount (1 −

_{ij}*c*)

*N*

^{−1}of synaptic weights is randomly picked and set to zero to represent sparseness. The connectivity scheme used here is in line with previous work outlining the importance of inhibitory processes in the generation of oscillatory cortical activity patterns at the mesoscopic scale (Wilson and Cowan, 1972; Amari, 1977), most notably with respect to alpha synchrony (Lopes da Silva et al., 1976; Mazaheri and Jensen, 2008, 2010; Lorincz et al., 2009) and slow-wave activity (Melzer et al., 2012).

Given our aim to assess oscillatory synchrony in the network, we were careful not to choose connectivity structures that would lead to localized states of persistent activity or propagating activity waves that are known to exist in systems such as Equation 1 (Roxin et al., 2005). We instead focus our attention on regimes for which the activity remains homogeneous (i.e., spatially uniform), which is ensured by the present connectivity scheme. We note that Equation 1 corresponds to a discretized neural field equation, for which the connectivity kernels would have been corrected for sparseness and randomness. In the absence of stimulation, i.e., ** I** = 0, and for strong and wide enough inhibition (i.e.,

*g*is large enough while

*r*is taken to be sufficiently small), the network in Equation 1 engages spontaneously in a stable alpha activity, oscillating with a frequency of ∼10 Hz, despite the random and sparse nature of the connections. The frequency of the oscillations is primarily determined by the mean delay τ, in a way we will detail below. Oscillatory dynamics of the network is shown in Figure 1

*C*.

#### Mean-field dynamics

The use of mean-field representations is a well established approach used to characterize the dynamics of large neural ensembles and is thus well suited to assess the impact of stimulation at scales that are relevant to population-level recordings in which spectral fluctuations might be observed, such as EEG (Deco et al., 2008). In the following, we develop a representation of driven network dynamics from the perspective of an ensemble average < ** u**(

*t*) > =

*ū*(

*t*) seen as a scalar readout measure of network collective activity when the population is large. Previous theoretical studies have shown that the ensemble average

*ū*(

*t*) is a reasonable model for the EEG (Robinson et al., 2001; Nunez and Srinivasan, 2006). The mean-field formulation further holds on presumed adiabatic dynamics, in which the system evolves in a mean-driven regime where the external fluctuations have smaller amplitudes compared with the autonomous dynamics of the system (i.e., the oscillations). The input timescale is also considered to be small compared with the synaptic rate α and interaction delay τ. The stochastic stimuli are also assumed to be sufficiently small and further to exhibit stationary statistics. Given that these conditions hold, the activity of individual units in Equation 1 can be expressed as small deviations

**of dimension**

*v**N*from the ensemble average

*ū*, i.e., where we introduce the unit vector 1. Since |

*v**| ≪ |*

_{j}*ū*|, we assume that the small deviations are subjected to the mean-corrected stimulus and obey the Langevin equation: where μ = <

**> stands for the ensemble average of the stimulus. Replacing Equation 2 into Equation 1 yields the following: Here, μ = <**

*I***> is the ensemble average of the stimulus. Deviations [**

*I***]**

*v**=*

_{j}*v*obeying the dynamics in Equation 3 above are identically distributed with probability density ρ

_{j}*(*

_{v}*v*) with zero mean and moments denoted <

*v*>.

^{n}Taking the ensemble average of Equation 4 yields a mean-field description of the dynamics of *ū*(*t*) when *N* goes to infinity. Given that *g* ∼ O(*N* − 1), we may use the independence and ergodicity of the local perturbations ** v** to obtain an expression of the corrected nonlinear neural response function (Shiino, 1987):
which can further be written as follows:
with the probability densities ρ

*(*

_{v}*v*) and ρ

*(*

_{w}*w*) and where <

*v*> =

^{n}*w̄*(

*c*) represents the mean effective connectivity and

*F*is the stimulus-corrected response function of the neurons.

The above rather technical derivations provide a clear message: the shape of the nonlinear response function of the neurons is altered by the inputs. As seen in Equation 5, the system response function now fully depends on the input statistics, and exhibits stimulus-dependent corrections to all nonlinear orders via the moments < *v ^{n}* >. Consequently, the external input alters the system's effective topology and dynamics.

##### Stimulus-corrected response function.

In the network, inputs to the neurons take the form of stochastic trains of actions potentials. We thus defined neural inputs as independent Poisson shot-noise signals of rate λ and amplitude *S*,
where the pulses δ(.) are Dirac delta functions with random arrival times *t _{i}^{k}* obeying a Poisson distribution. We assume that

*I*(

_{i}*t*) displays stationary statistics, i.e., both amplitude

*S*and rate λ are constant in time. We also assume that the fixed rate λ is large enough such that ρ

*(*

_{v}*v*) displays a stationary Gaussian profile. In contrast to other stochastic signals, such as Gaussian white noise, for instance, the mean and variance are not independent. It suffices to characterize the influence of the mean μ=<

**> =**

*I**S*λ on the dynamics, for which the variance of the linear process in Equation 3 becomes σ

_{v}

^{2}=

*S*

^{2}λ =

*S*μ (Gardiner, 2009). Using these expressions and inserting Equation 5 into Equation 4, the input-corrected mean-field dynamics simply reads For high values of the gain β, the response function

*f*behaves like a step function centered at zero, i.e.,

*f*(

*u*) ≈

*H*(

*u*) where

*H*is the Heaviside step function such that

*H*(

*u*) = 1 for

*u*≥ 0 and zero otherwise. In such conditions, the corrected response function

*F*in Equation 6 above simplifies to the following (to second-order): where erf [·] is the error function and σ

_{v}

^{2}= <

*v*

^{2}> =

*S*

^{2}λ =

*S*μ is the variance of the linear zero-mean process in Equation 3 with Poisson inputs. Once more, comparing Equation 6 and Equation 1 elucidates the influence exerted by the inputs on the recurrent architecture of the system. Figure 2

*B*shows the nonlinear gain, which is proportional to the system's response to external inputs. In sum, increasing the input variance renders the response function more flat and the stimulus response smaller.

#### Network stability

Based on the derived mean-field dynamics, the steady-state activity in Equation 6 is given by the following:
The equilibrium *ū _{o}* results from the interplay between the mean excitatory force delivered by the input, and the feedback it triggers in the network. As seen in Fig. 3

*A*, the equilibrium activity changes nonmonotonically as the Poisson stimulus intensity is increased. The stimulation train inhibits the network for moderate intensities, and then excites it as the intensity increases revealing how the input regulates mean network activity. Stimuli also alter the stability network dynamics. To see this, a linearization of Equation 6 about the equilibrium state given in Equation 8 yields the following: with the characteristic equation for the nonlinear susceptibility

*R*[

*ū*, σ

_{o}_{v}

^{2}] =

*g*[

_{o}w̄ F′*ū*], where

_{o}*F′*[

*ū*] is the nonlinear gain shown in Figure 2

_{o}*B*. The variable λ =

*a*+

*i*ω is the complex eigenvalue, where ω defines the frequency of the oscillations and

*a*is the damping rate. Figure 3

*B*shows that the susceptibility varies significantly whenever the input fluctuates and its amplitude |

*R*| first increases and then decreases as the stimulus intensity increases. This peculiar nonlinear mapping of the network response is determinant, as it indicates that the stimulus influences bidirectionally the efficacy of network interactions: changes in intensity can both increase or decrease the susceptibility.

This influence of stimuli on susceptibility also provides a mechanism underlying input-induced suppression of oscillatory activity, in which incoming inputs fully destabilize ongoing network oscillations, such as during event-related desynchronization observed in EEG. Indeed, using Equation 10, one can also determine the stimulus amplitude and/or rate at which stability of synchronous network states is lost. This occurs whenever the following equation
is satisfied for a particular set of stimulation parameters. A specific example of this can be seen in Figure 4*B*, where the stimulus causes a complete suppression of ongoing cyclic activity. Stimuli statistics in this simulation were chosen purposely such that they exceeded the stability threshold of endogenous network oscillations.

##### Driven oscillations.

Having established a mean-field representation of Equation 1 and exposed input-induced corrections to the system's nonlinear response function and stability equation, we now needed to expose the mapping between oscillation frequency, amplitude, and input statistics. Transition to synchrony in the network occurs in Equation 6 whenever Equation 10 above is satisfied for *a* = 0 and ω ≠ 0, in which case the system undergoes a so-called supercritical Andronov–Hopf bifurcation, transitioning from steady activity about *ū*(*t*) = 0 to global cyclic activity, or the other way around. The implicit dependence of the susceptibility *R* on input statistics, as seen in Figure 3*B*, implies that extrinsic fluctuations impact the oscillatory properties of the network. In more detail, an oscillatory instability (Andronov–Hopf bifurcation) emerges in the mean-field system for eigenvalues λ = *i* ω* _{c}* and
for critical values

*R*and delay τ

_{c}*.*

_{c}Aware of the effect of the inputs on the susceptibility *R*, perturbation analysis around the oscillatory bifurcation allows us to observe the influence of input statistics on the resulting peak frequency. Assuming that the system's baseline frequency is close to the critical frequency of the bifurcation, i.e., ω* _{o}* ≈ ω

*, effects of the stimuli can be seen as perturbations around the baseline, i.e., nondriven state. The peak frequency reads as follows: with γ = (2 ατ + 1/(ω*

_{c}_{c}

^{2}τ

^{2}− (1 + ατ)

^{2})), ΔR = (

*R*−

*R*)/

_{c}*R*. Equation 11 provides the desired mapping between input parameters and peak frequency and relates observed oscillatory activity with underlying stimulation dynamics. Fig. 4,

_{c}*C*and

*D*, displays the normalized frequency

*R*, cf. Fig. 3

*B*, whereas the amplitude saturates for large input stimuli.

In summary, the synergistic interaction between input dynamics and network nonlinear structure (which takes place via changes in the neural response function and in the susceptibility *R*) generically shapes the peak frequency by distorting the stability equation and provoking the observed change in oscillation frequency. Fluctuations in input mean and variance are echoed on-line by concomitant changes in susceptibility, causing the network peak frequency and amplitude to waver. The effect is also bidirectional: as plotted in Figure 3*C*, the spike trains delivered to the networks first slow down ongoing cyclic activity, and the normalized frequency displays a value smaller than 1. The trend changes, however, as μ = *S*λ increases and the peak frequency of the network starts to increase: stimuli now accelerate the dynamics. Together, these analyses show how spectral features can be modulated by stimuli without any change in circuitry.

## Results

To understand the contribution of input statistics on time-dependent fluctuations in power and peak frequency, we investigated the influence of irregular spiking inputs on the ongoing cyclic activity of a neural population using a combination of dynamical systems analysis and numerical simulations. Inputs to the network take the form of trains of spikes with time-variable amplitudes and rate, representing signals sent to the network by other neuronal populations located upstream, relaying changes in physical stimuli or alterations in cognitive state (Fig. 1*A*). As seen in Figure 4*A*, brief and strong stimuli generate stereotyped, transient responses predominantly in the alpha band. Ongoing oscillations remain stable after the stimulus and express identical peak frequency and power. However, when such pulses are prolonged as pictured in Figure 4*B*, suppression of activity can be observed. Yet, after stimulus offset, spectral features of the system are recovered. These more transient responses are to be contrasted with those generated by weaker, dynamically varying stimuli in which spectral modulation caused by inputs can be observed (Fig. 4*C*,*D*). Indeed, as seen in Figure 4*C*, changes in the stimulation intensity generate epochs of acceleration and deceleration, i.e., ongoing oscillations respond either by increasing or decreasing their frequency as the stimulus varies. Throughout, the amplitude of network oscillations fluctuates along with input intensity. This becomes even more transparent in Figure 4*D*, where the stimulus amplitude follows a sine wave. There, one can clearly observe the nonlinear mapping of stimulus response in the spectral domain. Yet, aside from the stimuli, no parameter was varied and this effect is purely input driven.

To understand the mechanism responsible for this effect, we developed and investigated a mean-field representation of the dynamics (see Materials and Methods), focusing on stimulation signals with stationary statistics, i.e., with nonfluctuating amplitude and rate. The stochastic analysis revealed that inputs support a gain control mechanism in which susceptibility is enhanced with respect to baseline in a nonmonotonic fashion. Based on our mathematical framework, it was possible to show that the nonlinearity of the network response, going from decelerated to accelerated synchrony, follows the stochastic dependence of the susceptibility: the frequency is proportional to |*R*| (cf. Fig. 3*C*; Eq. 11). The amplitude (and thus power) of network oscillations was also found to be tightly linked to the stimulus intensity, yet in a monotonic fashion. Figure 5 exposes the effect of the external stochastic input on the system's oscillation frequency in more detail. Increasing the amplitude of stimulation trains at weak intensities slows down ongoing activity by reducing the network oscillatory frequency (ω) with respect to the nondriven state, oscillating at baseline frequency (ω* _{o}*). The opposite occurs at larger intensities when the network peak frequency is shifted to higher values with respect to the nondriven state while increasing the stimulus intensity: the effect is thus bidirectional.

The analysis further exposed a nonlinear mapping between the stimulus intensity μ = Sλ, and the network oscillation frequency and its amplitude, pictured in Figure 6. The magnitude of the frequency modulation is substantial: it allows fluctuations of about 20% around ω* _{o}* induced by stimulus amplitude fluctuations of one order of magnitude. Strong input-induced decelerations were notably found to occur in a specific region of parameter space only, where stimulation rate is very high (100–500 Hz) but displayed amplitude is weak. We note that throughout, model parameters were kept constant and only the input was allowed to vary, showing the stimulus-driven nature of the effect.

## Discussion

Oscillatory activity is a key component of brain dynamics and has been the focus of an increasing number of neuroscientific investigations. For example, neuronal oscillations have been considered a possible mechanism through which internal states exercise top-down influences on stimulus processing to impact perception (Engel et al., 2001; Varela et al., 2001). In particular, the synchronization of oscillatory components seems to be relevant for many cognitive processes (Fell and Axmacher, 2011) and may operate across multiple scales of brain circuitry (Singer and Gray, 1995; Engel and Singer, 2001). In light of such, it is critical to understand how these oscillations are shaped by the inputs they receive and how this impacts information processing and communication across systems (Buzsáki and Wang, 2012). Time-dependent transitions in oscillatory activity are reliably observed across a wide range of frequencies in neurophysiologic measures (Lakatos et al., 2008; Schroeder and Lakatos, 2009; Ray and Maunsell, 2010; Thut et al., 2012; Van Zaen et al., 2013) and across a variety of computational models (Cohen, 2014). In particular, fluctuations in gamma peak frequency have been well characterized (Whittington et al., 1995; Wang and Buzsáki, 1996; Jadi and Sejnowski, 2014) while lower frequencies, which likely rely on different cellular mechanisms and operate over longer timescales, were also found to be significantly volatile (Haegens et al., 2014).

The robustness of these observations, along the short timescales on which they occur, suggest that peak frequency and power fluctuations may rely on a generic input-response feature of recurrent neural networks. The frequency and power/amplitude of coherent synchronous oscillations is commonly thought to reflect the underlying structural properties of neural networks and the features of its neural constituents. A direct consequence of this view is that the power spectrum of brain signals is oftentimes erroneously assumed to be time stationary. Our results add to a growing body of evidence showing that this is not the case, even in simplistic networks, and further add that the jittering of frequency and power must be regarded as a dynamic signature of fluctuating activation patterns. Using a fully nonlinear and stochastic approach, our analyses suggest that the spectral fluctuations observed in recurrent neural networks are caused by a synergetic interaction between inputs and the system's nonlinearity; an interaction that shapes the system's response function and resulting oscillatory behavior. While our study has primarily focused on the impact of afferent spiking inputs (by looking at Poisson shot noise), our conclusions would apply to a broader range of stochastic signals (e.g., Gaussian white noise) as well as those emanating from other, nonafferent sources. The range of frequencies expressed by our network, here located in the alpha band, can also be modified by adjusting the system's parameters, such that our findings would apply to other frequencies as well. We further add that despite the fact that our analysis was focused on input-induced spectral modulation, our model was nonetheless capable of expressing responses similar to those commonly observed with electroencephalography, such as event-related potentials (Fig. 4*A*) and event-related desynchronization (Fig. 4*B*).

Delay-induced synchrony is a candidate mechanism for the generation of oscillatory resting-state activity in the brain (Deco et al., 2011; Nakagawa et al., 2014), which is clearly distinct from the more local mechanism believed to be responsible for gamma activity. Operating at larger spatial scales for which synaptic conduction time lags become substantial (Cabral et al., 2014), our analysis shows that delay-induced synchrony can sustain the experimentally observed spectral variability of the alpha peak frequency, and thus that the power spectrum of neuroelectric signals can be seen as nonstationary across a wide range of frequency bands. As such, our work adds to a growing body of experimental and theoretical studies reporting complex, nonlinear interactions between endogenous oscillatory activity across multiple frequency bands and external stimulation patterns that build on and substantially extend various mechanisms such as resonance (Spiegler et al., 2011; Thut et al., 2012), plasticity (Fründ et al., 2009), or input-induced changes in the neurons' response functions (Doiron et al., 2001).

Cohen (2014) has recently put forward the concept of “frequency sliding” to describe dynamic changes in the peak frequency of modeled and empirical neural responses. The present complementary results go a step further by detailing a generic mechanism potentially responsible for these fluctuations and their reliability, providing a framework in which the spectral features of cortical networks can be directly linked to input statistics in a rigorous quantitative manner. In this vein are recent advances in adaptive frequency tracking, which have likewise provided empirical support to the proposition that stimuli can produce robust changes in the frequency of the responses of neural activity across a wide range of frequencies (cf. Van Zaen et al., 2013, their Fig. 8), which in turn vary as a function of whether or not stimuli induced the perception of an illusory contour. The directionality (increase/decrease) likewise varied as across frequency bands, which suggests that inputs might recruit differentially distinct recurrent structures. Such results provide additional empirical demonstration of the general principles revealed by the present computational model. Collectively, these show that network parameters alone cannot provide a full account for the observed pattern of neural responses; stimulus statistics and features are of equal relevance.

Response profiles of the modeled network followed a pattern consistent with divisive normalization (Carandini and Heeger, 2012), which is expected given the inclusion of diffuse, spatially wide inhibition in our model. Importantly, our results also suggest that mechanisms of divisive normalization are directly linked to those mechanisms engendering bidirectional frequency modulation. Such a link has also rightfully been outlined in both computational (Jadi and Sejnowski, 2014) and empirical (van Atteveldt et al., 2014) studies. Accordingly, our results show that normalization can be enhanced or suppressed, on-line, according to stimulus statistics and dynamics. This is also in accordance with the recent propositions of Cohen (2014) that frequency sliding may constitute a type of gain control mechanism wherein slower frequencies promote a lower threshold for firing while faster frequencies do the opposite.

We note that the present findings are distinct from resonance and entrainment phenomena in which rhythmic inputs either amplify (Ali et al., 2013) or override ongoing oscillatory activity effectively replacing the peak of the power spectrum (Thut et al., 2012). These mechanisms rely on the input itself having a rhythmic structure. Instead, our findings show how an unstructured stochastic input can bidirectionally modulate the power spectrum of ongoing oscillatory activity. Thus, in an era where neurostimulation becomes increasingly popular to support a wide variety of clinical interventions (Paulus, 2011; Baertschi, 2014; Cabrera et al., 2014; Gevensleben et al., 2014; Glannon, 2014; Schlaepfer and Bewernick, 2014), our results promote, in contrast, the idea that persistent yet weak stimulation of cortical networks can significantly amplify the power displayed at specific frequencies most likely engaged in the processing of sensory information, a mechanism that could constitute an alternative and/or complementary strategy to noninvasively tune ongoing cyclic brain dynamics.

### Conclusion

We show that the peak frequency and power displayed by neural populations evolving in coherent synchrony are context-dependent and nonstationary quantities that fluctuate according to the statistics of the inputs. This effect is due to a synergetic interaction between the system's intrinsic nonlinearity and input variability. As such, oscillatory properties of recurrently connected neural populations cannot be regarded as static, and the peak frequency must be seen as echoing the dynamics of the drive received by the neurons. The existence of these frequency transitions further suggests that classical neural spectral bands (e.g., theta, delta, alpha, beta, gamma, etc.) are mutually bound, where power can naturally transit from one band to another due to the action of exogenous and endogenous drive. New developments in adaptive frequency tracking are allowing for more accurate quantification of oscillatory activity and its dynamic modulation (Uldry et al., 2009; Van Zaen et al., 2010, 2013). Notably, spectral dynamics could be used alongside other methods to probe ongoing activity of cortical networks and trace event-related activation patterns in real time. Furthermore, our findings provide an important insight for future developments in quantitative analyses and simulations of oscillatory dynamics, including modeling of their alteration and breakdown to emulate disease states.

## Footnotes

This work has been supported by the Natural Sciences and Engineering Research Council of Canada (to J.L.), the Swiss National Science Foundation (Grant 320030-149982 to M.M.M. as well as the National Centre of Competence in Research project “SYNAPSY, The Synaptic Bases of Mental Disease” [project 51AU40-125759]), the Swiss Brain League (2014 Research Prize to M.M.M.), and the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 257253 to A.H.).

The authors declare no competing financial interests.

- Correspondence should be addressed to Dr. Jérémie Lefebvre, 53 rue Sainte-Marie, Gatineau, Québec, Canada, J8Y2A6. jeremie.lefebvre{at}hotmail.com