Poisson-distributed noise induces cortical γ-activity: explanation of γ-enhancement by anaesthetics ketamine and propofol

Additive noise is known to affect the stability of nonlinear systems. To understand better the role of additive noise in neural systems, we investigate the impact of additive noise on a random neural network of excitatory and inhibitory neurons. Here we hypothesize that the noise originates from the ascending reticular activating system. Coherence resonance in the γ-frequency range emerges for intermediate noise levels while the network exhibits non-coherent activity at low and high noise levels. The analytical study of a corresponding mean-field model system explains the resonance effect by a noise-induced phase transition via a saddle-node bifurcation. An analytical study of the linear mean-field systems response to additive noise reveals that the coherent state exhibits a quasi-cycle in the γ-frequency range whose spectral properties are tuned by the additive noise. To illustrate the importance of the work, we show that the quasi-cycle explains γ-enhancement under impact of the anaesthetics ketamine and propofol as a destabilizing effect of the coherent state.


Introduction
Complex systems exhibit a hierarchy of scales, possibly both in space and time.Microscopic system elements interact and may synchronize and thus represent a functional (abstract) entity on a larger mesoscopic scale.Several of these mesoscopic functional entities may interact as well and build abstract macroscopic functional units.This hierarchy [1] may exhibit a large number of levels.Examples of such scales can be found in biology as animals moving in swarms, in sociology as humans aligning their behavior to a certain trend or in physics as fluid molecules moving in a coherent structure.
There is increasing evidence that CR plays an important role in brain [12,19].In this context, the question arises where the noise originates and how it is generated.Besides intrinsic fluctuations on several scales, such as ion channel fluctuations and synaptic noise, recent work has hypothesized that the impact of the ascending reticular activating system (ARAS) [20] may be modeled as a noise source [21][22][23].This hypothesis describes the activity of the rather complex ARAS network of brain areas by unspecific random activity.The hypothesis states that enhanced ARAS activity has a similar impact as an increased noise level.Since the ARAS is known to set the level of excitation of brain areas, high ARAS activity can be modeled by a high input noise level.A previous work [24] has provided some experimental support for this hypothesis demonstrating that noisy neurostimulation may enhance the level of brain excitation.Moreover, there is some experimental evidence that the ARAS tunes the mammalian attention and affects gamma-activity in visual cortex [25,26].As a consequence of these pieces of evidence, we hypothesise that ARAS tunes the brain by CR.
To this end, the current manuscript presents a study of a random neural network subjected to additive uncorrelated noise.It extends a previous study [27] by several elements.We show CR in the γ-frequency range (30-60 Hz) numerically for a broad range of Poisson noise levels and explain the underlying mechanism analytically by a mean-field model.We find a new second bifurcation for very large noise levels that complements previous known CR results nicely.An additional detailed analytical linear response study reveals the noise impact on the system's coherent state.We identify a quasi-cycle as underlying process of the coherent state, provide a short analytical description of these quasi-cycles and show how noise tunes the spectral power of the coherent state.In addition, the γ-activity enhancement found experimentally in electroencephalographic (EEG) under anaesthetic medication is explained by a Hopf bifurcation of the noise-induced coherent state.

Material and methods
At first, we present the network model and then show the derivation of a mean-field model that captures the major dynamics of the network.A short analytical study of the systems topology complements this part.In a subsequent section, a short response theory of the system close to a single equilibrium describes the linear evolution about the equilibrium and the systems power spectrum.

The network model
We assume a network of excitatory and inhibitory neurons with N neurons in each network.Excitatory (inhibitory) neurons excite (inhibit) each other according to the connectivity matrix F ∈ R N×N , and they excite (inhibit) inhibitory (excitatory) neurons according to the connectivity matrix M ∈ R N×N .These neural network interactions are described by a dynamical system governing the evolution of the state variable vectors Such a model is a typical rate-based models [28].It implies asynchronous activity and synaptic response functions of first order with time scales τ e,i .The state variables V and W represent excitatory and inhibitory dendritic currents, respectively, while ξ e,i ∈ R N are stochastic inputs from various sources, such as ion channel fluctuations or stochastic input from other brain areas such as the ARAS.Moreover, I e,i is constant input with e = (1, . . ., 1) t and S 1,2 [u] ∈ R N is the nonlinear transfer function with (S 1 [u]) n = H 0 H(u n ), (S 2 [u]) n = H(u n ), H 0 > 0 and the scalar Heaviside transfer function H with H(u) = 0 ∀u < 0, H(u) = 1 ∀u 0.
The interaction networks are directed non-sparse Erdös-Rényi networks with connection probability density c = 0.95.It represents a reasonable network topology on a microscopic scale, i.e. in patches of diameter below 1 mm [29].We assume non-vanishing matrix elements as (F) ij = F 0 /cN and (M) ij = M 0 /cN.Then the edge spectrum of F contains its maximum eigenvalue and eigenvector [27,[30][31][32][33] and its bulk spectrum has the maximum eigenvalue Equivalent results hold for matrix M. It is obvious that λ 2 λ 1 and λ 2 ≈ 0 for large mean degree cN.If N or c decreases, then λ 2 increases and the spectral gap becomes smaller and this approximation does not hold anymore.As a consequence, the system's stability may change, cf appendix of [34].
We assume the noise process at excitatory neurons {ξ e n (t) } to have the mean ξe and variances D 1 with The input to a neural population is well-described by incoming spike trains that induce dendritic currents at synaptic receptors.In a good approximation, the neurons interspike interval obeys a Poisson distribution [35] and incoming spike trains at mean spike rate λ in induce random responses at excitatory synapses with time constant τ in .This random process has the ensemble mean and ensemble variance [36,37] in λ in τ in /2 assuming the synaptic coupling weight w in .Since a Poisson distribution converges to a Gaussian distribution for large enough mean, we implement this input current as a Gaussian random process with the corresponding mean and variance while ensuring the validity of this approximation by a large enough input firing rate λ in .It is important to point out that both mean and variance are proportional to the input firing rate.
Moreover, the noise process at inhibitory neurons (ξ i ) n at network node n is Gaussian distributed with zero mean, noise intensity D 2 and uncorrelated in time Parameters chosen for the network simulations are given in table 1.

The mean-field model
In the model (1), the system activity V is expanded into the biorthogonal basis with modes Then the system activity is decomposed by a n Φ e n with complex mode amplitude a n ∈ C. Here, † denotes the transpose complex conjugate.The same decomposition is done for W with the basis with the complex mode amplitude b n ∈ C and the biorthogonal basis Now the projection of the systems activity V, W onto the respective basis {Ψ e k } and {Ψ i k } yields the amplitude equations Choosing Ψ e k , Φ e k as eigenvectors of F with eigenvalue λ e k ∈ C and and In addition ) with m e,i (t) = e t ξ e,i (t)/N.Equations ( 6) and ( 7) describe a multivariate Ornstein-Uhlenbeck process whose solutions obey for t → ∞.Moreover the terms V, W in equations ( 4) and ( 5) can be written as [27] V(t) = (a 1 (t) + s e (t))e + w e (t) with e −(t−τ )/τ e,i η e,i (τ )dτ ( 10) e −(t−τ )/τ e,i ρ e,i (τ )dτ (11) with η e,i (t) = ξ e,i (t) − e ξe,i 0 , e t η e,i (t) = Nρ e,i (t) with ρ e,i ∼ N (0, D 1,2 /N) and temporally constants ξe,i 0 .The term ρ e,i represent finite size-fluctuations with variance D 1,2 /N that vanish for infinite size networks, i.e. ρ e,i → 0 for N → ∞.
The mean-field equations can be written as + I e + ξe 0 + ρ e (t) Since w e,i in equation ( 10) is a multivariate stationary Ornstein-Uhlenbeck process, it is ergodic and its statistics over time is identical to the ensemble statistics at a single time point.Hence w e,i is stationary over the network with a stationary probability density function Then where the approximation is good for large N. Similarly, Essentially, for the Poisson noise s e (t) = ρe (t)e with finite-size fluctuations ρe (t) ∼ N (0, D 1 /N 1 ), ξe 0 = w in r in τ in and D 1 = w in ξe 0 /2 yielding the mean-field equation We observe that the mean-field equation is driven by a stochastic multiplicative force ρe,i (t) and an additive force ρ e,i (t).These stochastic forces result from the finite network size and vanish for N → ∞.
To understand better the dynamics of the system, we compute the equilibrium (a 1 0 , b 1 0 ) for vanishing stochastic forces and analyze its noise-free stability.The equilibrium is the solution of Please note that the equilibrium does depend on statistical moments of the additive noise.By virtue of the connectivity kernel choice, they are independent from the network size N and the connectivity probability c.Then linearizing the system (15) about the equilibrium with small deviations x(t) = a 1 (t) − a 1 0 and with the constants G 1 = dG 1 da 1 and G 2 = dG 2 db 1 computed at the corresponding steady state.Here, we neglect the multiplicative stochastic force.This differential equation system can be written as with The eigenvalues λ of the matrix A are defined by Av = λv with eigenvector v and are found as They provide information on the stability of the noise-less system around the chosen equilibrium.If the trace Tr(A) < 0, then the equilibrium is asymptotically stable.

Linear response: power spectrum
The following paragraphs assume stability of the equilibrium and the system evolves in a stationary state.We consider the linear system (18) driven by uncorrelated noise ξ(t).In the stationary state, we can express the inhomogeneous solution of equation ( 18) for t → ∞ as where G ∈ R 2×2 is the Green's function.Its Fourier transform G(ω) is defined as with the identity matrix I ∈ R 2×2 .Since the system is in a stationary state, the Wiener-Khinchin theorem applies and the spectral decomposition of a solution component z(t) = X i (t), i = 1, 2 is given by its power spectrum where Here and in the following, * denotes the complex complex-conjugate and • represents the ensemble average.Equation (22) holds true since the system is stationary in time and thus is ergodic.Then rewriting equation (20) as and inserting this equation in equation (22) we gain

we can then rewrite the correlation function as
and, utilizing equation ( 21), gain the power spectrum More specifically, the spectral power density of X 1 (t) = x(t) in equation ( 17) reads with the matrix elements A 11 , A 12 of matrix A. Please re-call, that Di ∼ 1/N and hence S 2 (ω) → 0 for N → ∞.
In other words, the system does not evolve for infinitely large networks.Conversely, the smaller the network the larger the macroscopic finite-size noise.

Network activity
Numerical simulations of equation (1) reveal CR. Figure 1(A) show a jump from a non-synchronous noisy state at low Poisson input rate λ in to a coherent state at larger input rate.Figure 2 shows the corresponding power spectrum in each state with a clear spectral peak in the coherent state.Further increasing the Poisson noise input rate (figure 1(B)), coherent oscillations become less prominent and the system exhibits uncorrelated activity again for large input rates.
In addition to the coherence, increasing the input rate the system transitions from a high activity state over a low-activity state (that is coherent after the transition) at intermediate input rate back to a high-activity state for large input rates.
The networks mean dynamics does not depend on the network size.Figure 1(A) demonstrates that the critical input rate, at which the transition to a coherent state emerges, does not depend on the network size.However, the smaller network (N = 50) exhibit stronger finite size fluctuations around the mean, than the large network (N = 200).

Mean-field: bifurcation analysis
The findings in the network may be explained by the mean-field description (15).To gain a first insight, we have computed the equilibria and their stability of the system.The equilibria are solutions of equation (16).Figure 3(A) shows a 1 0 subjected to the Poisson input rate λ in .For small input rate, the topology of the system exhibits a stable upper state, a saddle center state and a lower state that is either a stable or unstable focus.This explains the high-activity state observed in figure 1.For larger input rates, the system evolves about a stable focus equilibrium in a lower state.This state is observed in figure 1 as well and represents the coherent low-activity state.Hence, the system can be bistable, meaning that both the upper and the lower states are present and stable.In this case the neurons activity in the stochastic system can randomly switch from one state to the other.This switching may occur for large noise input and are observed in figure 1(A) for N = 50.There, the macroscopic noise fluctuations ρ e,i (t) in the mean-field are large finite-size fluctuations and may induce these jumps in the bistable state.However, here we have to point out that for such small networks the validity limits of the mean-field description may be reached and this provided explanation may be valid qualitatively only.The system also may be monostable, i.e.only the upper (for small and large noise input rate) or the lower state (for intermediate noise input rate) is present and stable.Taking a closer look at the low-activity state, we observe a Hopf bifurcation with critical frequency in the γ-frequency range (see figure 3(B)).Increasing the noise input rate the equilibrium bifurcates from an unstable focus to a stable focus and its eigenfrequency diminishes down to zero when the state ceases to exist.The subsequent section elucidates further the spectral properties of this lower state.For low input rate, the system is monostable at an upper equilibrium branch, then close to the bifurcation point, the system becomes bistable with a stable node and a stable focus.For high input rate, after the bifurcation point, the system becomes again monostable and is a stable focus.(B) Maximum real part of the eigenvalue (black) and the eigenfrequency (red) subjected to the input rate λ in in the lower equilibrium branch.We observe a Hopf bifurcation at ∼1000 Hz.This topology describes well the CR: increasing the noise level with an initial state about the upper equilibrium, the system performs a noise-induced phase transition from the upper stable node to the lower stable focus via a saddle-node bifurcation.Then further increasing the noise there is a second saddle-node bifurcation and a noise-induced phase transition from the lower stable focus to the upper stable node.

Mean-field: the linear response of the coherent state 3.3.1. Quasi-cycles
To understand the dynamics of the phase-coherent state in the lower bifurcation branch in figure 1, we investigate the linear response of the system to the macroscopic finite size-fluctuations.Here, we point out that the mean-field itself is determined by the microscopic input noise ξ e , ξ i in equation ( 1), whereas the macroscopic fluctuations ρ e , ρ i echo this additive microscopic noise and remain as additive macroscopic finite-size fluctuations.Since these finite-size fluctuations represent the spatial average of the input noise and the full spatial set of microscopic noise processes tune the mean-field, we consider the macroscopic noise being independent of the microscopic noise.
The coherent state observed in figure 1 evolves about an equilibrium (a 1 0 , b 1 0 ), that is a stable focus if the eigenvalues of A, cf equation (19), have non-zero imaginary parts and its real parts are negative, i.e.
Then the corresponding angle frequency is Let us rewrite equation (25) as with the susceptibility function [38] defined as If the observed variable x is driven by weak noise, i.e.D1 ≈ 0, then the shape of the power spectrum is primarily determined by the susceptibility function χ.Then, the frequency of the spectral peak is defined by and the spectral density at the peak is At first, we observe a spectral peak at angular frequency ω s = 0, i.e. the system oscillates with a primary angular frequency ω s although the corresponding deterministic system is stable and does not exhibit an oscillation for t → ∞.The oscillation observed is induced by the additive noise and is called a quasi-cycle in the literature [39].In (29), f2 depends on D1 and D2 , which themselves are inversely proportional to N, rendering the spectral density S 2 also inversely proportional to N. For weak noise, the frequency of the quasi-cycle is close to the eigenfrequency ω l of the corresponding linear system, i.e.
Moreover, the smaller the eigenvalue real part Tr(A) in equation ( 29) the larger is the peak spectral power.In other words, an increase of spectral power of the peak reflects convergence to the system's stability threshold.This effect is essential to understand the drug impact discussed in the subsequent section.

Noise tunes γ-power
To illustrate how the additive noise tunes the power of the coherent state, figure 4 shows the spectral power distribution for a large range of input rate λ in .From an initial state at large noise input rate λ in in the lower equilibrium branch, decreasing the input rate increases the slightly the frequency of the spectral peak in the γ-frequency band and the corresponding peak power, cf figure 3(B).At ∼800 Hz, the system jumps from the lower equilibrium branch to the upper branch (cf figure 3), where the system does not oscillate anymore.
The power surge observed when decreasing the noise level can be understood by equation (29): the closer the system to the stability threshold, the smaller the eigenvalue real-part Tr(A) 2 /2 and the larger the peak power.This is consistent with figure 3(B).

Effects of ketamine and propofol
Modulation of cortical γ-activity in the brain may reflect modulation of attention or arousal.For instance, in general anaesthesia some drugs induce an enhancement of γ-activity in experimental mammalian EEG data, e.g. after the administration of ketamine [40,41] and propofol [42].Since our model describe interacting neural networks, it is interesting to investigate whether it may explain this γ-enhancement.(25).Both ketamine ((2) in figure 5) and propofol ((3) in figure 5), respectively color-coded in blue and red, cause an increase of the amplitude of the gamma peak, compared to the initial state ((1) in figure 5) color-coded in black.
To this end, we consider the experimental finding that propofol increases the decay time constant τ i of GABA A synapses [43], associated, in our model, with the inhibitory neurons population.For our model we will then simply study the action of propofol by studying the effect caused by a relative increase of τ i .For ketamine, it has been argued that it causes a decrease in the decay time constant of glutamatergic neurons [44], which we model here as an decrease of τ e .However, the effects of ketamine rely on a rather complex physiological mechanism, which might also cause an increase of τ i .In the following we will consider pure τ e effects and associate it with ketamine.
A detailed linear stability study of the lower equilibrium dependent on the synaptic time scale τ e and τ i reveals a Hopf bifurcation (figure 5) for which Tr(A) = 0.
Both decreasing τ e and increasing τ i , the system approaches the stability threshold and the γ-power increases, as seen in figure 6.This corresponds well to experimental findings [40][41][42].

Discussion
The present work has demonstrated CR in the γ-frequency range and explains the underlying mechanism analytically by a mean-field model.At low input rates, we observe a saddle node bifurcation at the upper branch and a Hopf bifurcation at the lower branch.The bifurcation thresholds are different in both branches which reflects a hysteresis.This can be observed by comparing figures 1 and 4: the former presents a different critical input rate at the upper branch than at the bottom branch shown in the latter figure.Moreover, we find a new second bifurcation for very large noise levels that provides a full picture of the coherent resonance phenomenon: for low and high noise levels the system exhibits non-coherent behavior while intermediate noise levels render the system more coherent.An additional detailed analytical linear response study reveals the noise impact on the system's coherent state and we identify a quasi-cycle as underlying process of the coherent state.Moreover, we demonstrate that noise tunes the power and frequency of the γ-activity in the coherent state.For evaluation, we show that the proposed noise-driven cortical model permits to explain the γ-power modulation in the presence of the drugs ketamine and propofol as known from literature.
The model under study simplifies interactions between neurons by a rate-coding assumption and assumes simple threshold firing rate models.Such a model does not capture diverse dynamical features of single neuron activity, such as voltage-dependent synaptic response, synaptic adaption or different firing types of neurons, e.g.neurons type I and type II.Previous work has involved such single neuron properties in rate-coding models [36,45,46].Future work will replace the Heaviside firing rate function in equation ( 1) by a more realistic model.This extension will not constrain the provided analysis.
Moreover, the gained results just indicate that the model will be applicable to explain several γactivity phenomena, but only the effect of ketamine and propofol has been shown.Future work will study other phenomena like γ-bursts [47] or the relation between attention, ARAS activity and γ-activity modulation [25].
Now we specify ξ(t )ξ T (t ) = 2δ(t − t )D D = D1 0 0 D2 , where D1 = D 1 /τ 2 e N, D2 = D 2 /τ 2 i N is the noise variance of the ith component of ξ.Inserting these expressions back into the correlation function and utilizing the Fourier transform of the Green function g(ω) yields

Figure 1 .
Figure 1.Increasing additive noise induces phase transitions from a non-oscillatory non-coherent state to an oscillatory coherent state (CR).Numerical simulation result of equation (1).(A) Phase transition to a coherent state for two network sizes.(B) First phase transition from a non-coherent to a coherent state at low input rates and at large input rates from a coherent state to a non-coherent state.

Figure 2 .
Figure 2. Non-coherent states exhibit a flat spectral distribution, whereas the coherent state shows a clear spectral peak in the γ-frequency range.The spectral distributions are computed from network average time series at λ in = 700 Hz (black), λ in = 1900 Hz (red) and, λ in = 9000 Hz (blue).The power spectra were computed by a Fourier transform from simulated time series of duration 5 s.

Figure 3 .
Figure 3. Bifurcation diagram of the mean-field equilibria.(A)For low input rate, the system is monostable at an upper equilibrium branch, then close to the bifurcation point, the system becomes bistable with a stable node and a stable focus.For high input rate, after the bifurcation point, the system becomes again monostable and is a stable focus.(B) Maximum real part of the eigenvalue (black) and the eigenfrequency (red) subjected to the input rate λ in in the lower equilibrium branch.We observe a Hopf bifurcation at ∼1000 Hz.

Figure 4 .
Figure 4. Noise tunes frequency and instantaneous power of γ-activity.The power spectral density is shown on a logarithmic scale.The network has been simulated numerically for 55 s with time step Δt = 0.05 ms while decreasing λ in linearly from maximum to the minimum rate inputs starting from the lower equilibrium branch.The spectral power distribution has been computed by a sliding Fourier transform with window size 2 s and shift time 0.5 ms.

Figure 5 .
Figure 5. Hopf bifurcation in the lower equilibrium dependent on synaptic scales τ e and τ i .The dashed line represents the stability threshold.The parameter space point (1) in black corresponds to the initial position in the τ e,i space.The effect of ketamine modeled by a decrease of τ e corresponds to a transition from point (1) to point (2) in blue.The effect of propofol modeled by an increase of τ i corresponds to a transition from region 1 to region 3 in red.The state approaches the stability threshold in the direction of decreasing τ e and increasing τ i yielding a power surge in the γ-range.

Figure 6 .
Figure 6.Ketamine and propofol both increase the amplitude of γ-oscillations in the analytical model.Ketamine and propofol effects here are modeled respectively by a decrease of τ e and an increase of τ i .The plotted functions represent the linear analytical power spectra(25).Both ketamine ((2) in figure5) and propofol ((3) in figure5), respectively color-coded in blue and red, cause an increase of the amplitude of the gamma peak, compared to the initial state ((1) in figure5) color-coded in black.