Phase Coherence Induced by Additive Gaussian and Non-gaussian Noise in Excitable Networks With Application to Burst Suppression-Like Brain Signals

It is well-known that additive noise affects the stability of non-linear systems. Using a network composed of two interacting populations, detailed stochastic and non-linear analysis demonstrates that increasing the intensity of iid additive noise induces a phase transition from a spectrally broad-band state to a phase-coherent oscillatory state. Corresponding coherence resonance-like system behavior is described analytically for iid noise as well. Stochastic transitions and coherence resonance-like behavior were also found to occur for non-iid additive noise induced by increased heterogeneity, corresponding analytical results complement the analysis. Finally, the results are applied to burst suppression-like patterns observed in electroencephalographic data under anesthesia, providing strong evidence that these patterns reflect jumps between random and phase-coherent neural states induced by varying external additive noise levels.


INTRODUCTION
Noise is omnipresent in nature [1] playing a synergetic role in numerous organic and inorganic processes [2]. Random spatial and/or temporal fluctuations in neural systems are essential for brain function: noise is a key component of numerous neural processes, allowing brain networks to enrich their dynamical repertoire [3], enhancing sensory processing and detection in animals [4], improving signal-to-noise ratio [5], and even by sustaining cognitive brain functions [6].
Much insight about how noise shapes neural processes can be gained by studying how noise impacts non-linear dynamical systems, especially their oscillatory properties. In contrast to multiplicative noise, which reflect parametric fluctuations that are known to transform the stability of numerous systems [7,8], additive noise describes fluctuations originating from unknown sources or external fluctuations from outside of the system. Previous studies of large-scale systems have provided multiple examples of additive noise-induced phase transitions [9][10][11][12]. Recent studies on the impact of additive noise on networks of simplified delayed non-linear oscillators have demonstrated that noise shapes non-linear network interactions, hence tuning its spectral features [13,14]. However, most of these studies have focused on dynamics close to instability [15] or on dynamics in weakly coupled networks [16]. It remains poorly understood how strong noise impacts natural systems that are strongly connected and far from the stability threshold [8,17,18].
To extend these results and explore how strong additive noise impacts a more realistic network and its oscillatory properties, the present work considers a generic non-linear activation-inhibition network model built of excitable systems. In addition to neural systems [19], such networks have also been observed in enzymatic systems [20] and autocatalytic patterns in sand dunes and other inorganic patterns [21]. In contrast to previous studies, here the isolated excitable populations do not oscillate on their own, and oscillatory behavior results from collective interactions tuned by additive noise. We show how additive (independent identically distributed or iid) Gaussian and (independent non-identically distributed or non-iid) non-Gaussian network noise shape non-linear interactions in the model, specifically how it impacts its oscillatory dynamics. Both iid and non-iid noise induce transitions from noncoherent non-rhythmic states to coherent oscillating states. To describe the observed phase transitions, we derive analytically a deterministic order parameter equation [17] where additive noise transforms the dynamics through a novel noise-dependent nonlinear response function. Using these results, we demonstrate and fully characterize coherence resonance-like behavior for varying iid noise and varying heterogeneity of non-iid noise. In order to gain the analytical results, we assume homogeneity of the networks, ergodicity in groups of neurons and statistical independence of connectivity weights and system dynamics. Taken together, our work provides convincing support for burst suppression [22] dynamics observed in the mammalian brain under anesthesia gathered through electroencephalography. Indeed, our model and analysis shows that a decrease in additive noise triggers a phase transition whose spectral signature provides a physiological support of phenomena seen experimentally.

METHODS
The subsequent sections show the derivation of mean-field dynamics and random fluctuations about the mean fields. The mean-field dynamics is determined by additive noise properties.

The Model
The present work considers a system of two networks with nonlinear activation-inhibition interactions. The networks under study exhibit random connections motivated by network models for associative memory in the brain [23,24].
We consider two coupled homogeneous networks V and W of number N whose nodes obeŷ for n = 1, . . . , N, V n ∈ V, W n ∈ W and the temporal operatorL = d/dt +Î, whereÎ is the unity operator. This model resembles the basic microscopic network model that yields the famous population model of Wilson and Cowan [25] and that can show oscillatory dynamics [26]. The network connectivity matrices F ∈ ℜ N×N and M ∈ ℜ N×N are random with matrix elements F nm , M nm ≥ 0. Since the network V excites other cells and W inhibits other cells, the network V and W is excitatory and inhibitory, respectively. Each network is homogeneous with for all n ∈ [1; N] and constants F 0 , M 0 > 0. This is a realistic assumption that has been affirmed experimentally [27]. In the following examples, without contraining generality we choose Erdös-Rényi networks [28] for simplicity with connection probability c, mean degree Nc, and connection weights

Additive Noise
Both populations in Equation (1) are subjected to additive noise ξ n (1,2) . However, we allow one population to be driven by iid noise (W), while the other (V) can be driven by either iid-or noniid noise. Specifically, according to this, and given that the noise processes are independent in time and between nodes we have for all n, m = 1, . . . , N with mean valuesξ (1) n , noise intensities D (1) n , D (2) and · denotes the ensemble average. We point out that all stochastic processes in network W have zero mean and identical variance D (2) and thus are independent and identically distributed (iid). In contrast, additive noise in network V may be inhomogeneous over nodes with node-dependent mean and variance. This case involves a superposition of different Gaussian distributions in the network and the network noise is nonidentically distributed and thus non-iid.
Specifically, each noise process in network V belongs to a certain class G m , m = 1, . . . , M of M classes. Noise processes ξ (1) n in a class G m , i.e., n ∈ G m , share their meansξ m and variances D m while different classes may exhibit different means and variances.
To derive the probability density function p V (ξ ) of the noise process over the network V at a certain time instance we assume that the noise processes in each class G m are ergodic, i.e., all noise processes in each class are stationary and hence the ensemble average equals the average over nodes. Then the mean and variance of class G m is where N m is the number of nodes in class G m . For the network mean with the constant weights p m = N m /N. Similarlȳ utilizing (5) and (6). These expressions for network mean and network variance show sums over the mean and variances of the noise classes {G m }. This indicates that the probability function of the network noise p V (ξ ) is a function of the probability density functions of the single classes as well. Indeed, we find For completeness, we define the average of additive noise in network W to N n = 1 ξ (2) n (t)/N =ε 2 .

Mean-Field Dynamics
By virtue of the random network connectivity, no finite spatial scale is prominent and the system converges to the mean field solution. In addition the noiseless case D (1,2) n = 0 stipulates a single constant stationary solution V 0 Consequently, we expect that a solution, that is constant over the network, describes the network dynamics. A mathematical proof of this assumption would exceed the current work and we refer the reader to forthcoming work. However, as will be seen in later sections, this assumption holds true in the network under study and permits to describe the network dynamics very well.
We choose the ansatz whereV,W are the averages over the respective network activity and υ n , π n denote the deviations from the average activity, with the mean deviationsῡ,π . Taking the sum N n = 1 /N on both sides of Equation (1) we find where e 1 (t) =ε(t) −ῡ(t) − dῡ(t)/dt, e 2 (t) =ε 2 (t) −π(t) − dπ(t)/dt. We point out thatε,ε 2 ,ῡ, andπ do vanish for N → ∞ due the statistical law of large numbers and do not vanish for finite-size networks. These values vary randomly in time and hence e 1,2 (t) are random as well. In the following, we assume large networks with |e 1 (t)| ≪ |I (1) 0 |, |e 2 (t)| ≪ |I (2) 0 |, and hence e 1,2 are considered as small stochastic forces that do not affect the dynamics of the mean-fields.
Since F nm and υ m are independent random numbers Frontiers in Applied Mathematics and Statistics | www.frontiersin.org utilizing the statistical independence between connectivity and stochastic process with the expectation value E m [X m ] ≈ N m = 1 X m /N. We point out that condition (12) holds for a large enough number of connected nodes, i.e., for large enough mean degree cN ruling out sparse networks.
Then condition (12) leads to and with the probability density function p V (υ, t) of the deviations υ at time t. The term ε represents the finite-size error of the integral. Equivalently with the estimation error η and the probability density function p W (π, t) of deviations π. This yields the evolution equations for the average activitŷ with F = N n = 1 F 0 /N = F 0 , M = N n = 1 M 0 /N = M 0 , and I 1,2 = I 1,2 0 + e 1,2 .
Here, it is important to mention that the deviations υ n , π n decouple from the averagesV,W and Equation (17) (11) and (18). If this does not hold true, υ n and π n depend onV,W and Equation (17) are not sufficient to describe the system dynamics.
The evolution equations (17) depend on the probability density functions p V (υ, t), p W (π, t). For large times, the solutions of (19) are [29] Utilizing the noise properties (3) we gain whereξ m , D m ,ε andD are taken from Equations (5)- (7). Then the probability density function of network V is Additive noise in network W is iid and hence p W (π) = N (0, D (2) ).

Power Spectrum
Power spectra of the average activities provide evidence on phase coherence in the interacting networks. To this end, we have computed numerically the power spectra of simulated activitȳ V(t) by the Bartlett-Welch method with time window 100 s and overlap rate 0.995. This numerical power spectrum involves the full non-linear dynamics of the system. To compute analytically the power spectrum of the average activity, we can approximate the nonlinear dynamics by considering the stochastic system evolution about an equilibrium. To this end, we can write Equation (17) as with the state vector Z = (V,W), the non-linear vector function N ∈ R 2 and the noise term e = (e 1 , e 2 ). Neglecting the noise, the equilibrium Z 0 = (V 0 ,W 0 ) is defined bȳ Then small deviations from the equilibrium u = Z − Z 0 obeys du/dt = Lu + e with the linear matrix and the linear power spectrum ofV(t) reads [30] R(ν) = D 0 (L 2 22 + L 2 12 + 4π 2 ν 2 ) 4π 2 (L 11 + L 22 ) 2 ν 2 + (detL − 4π 2 ν 2 ) 2 (23) with frequency ν, the determinant detL and the matrix elements L mn of L in (22). The terms S ′ 1,2 (x) = dS 1,2 /dx are the non-linear gains computed at the corresponding equilibrium. Since the noise e(t) results primarily from finite-size effects, its properties are hard to model. In Equation (23), we have assumed e 1,2 = 0, e k (t)e l (τ ) = D 0 δ kl δ(t − τ ) for simplicity. Consequently, the linear analytical power spectrum (23) does not include all nonlinear stochastic effects and does not capture the correct noise strength, but it conveys a good approximation of the numerical power spectra while involving the additive noise-induced shaping of the non-linear interactions.

Phase Coherence and Bandpass Filter
The Phase Locking Value (PLV) has been developed by Lachaux et al. [31] and estimates the phase coherence between two signals s 1,2 (t). The major idea is to compute the time-dependent phase, e.g., by wavelet transform, of each signal φ 1 (f , t) and φ 2 (f , t) in the frequency band about center frequency f at time t. This approach follows the idea of an instantaneous amplitude and phase of a signal (see e.g., [32][33][34][35] for more details). The implementation is based on the general idea of a linear filter where a signal s(t) is filtered by a sliding window function w. This window exhibits an amplitude-modulated complex-valued oscillatory function with intrinsic frequency f . Heuristically, the phase of the filtered signals at time t and frequency f is the mean phase of the signal at time t with respect to the phase of the window function w. Then one computes the circular variance (or PLV) in a time window of width w about each time instance t n . If R 1,2 (f , t) = 1 both signals are maximum coherent in the frequency band about f at time t, R 1,2 (f , t) = 0 indicates no phase coherence. This is the classical bi-signal phase coherence. To obtain a measure for global phase coherence as considered in the present work, one may average the PLV of two signals s n , s m over a set of phase pairs M of mass M [35] and gain the global Phase Locking Value that represents the degree of phase coherence in the system. Similar to the two-phase PLV, gPLV(f , t) = 1 reflects maximum phase coherence, i.e., all nodes in the system have identical phases, and gPLV(f , t) = 0 indicates that all nodes are out of phase. In this context, it is important to note that instantaneous phase values are defined properly only in a certain time-frequency uncertainty interval for non-vanishing instantaneous amplitude in the same interval. To avoid spurious phase coherent values, we consider those phase pairs only, whose instantaneous spectral power exceeds 50% of the maximum spectral power at the corresponding time instance [36].
Typical experimental electroencephalographic data (EEG) is bandpass filtered to remove artifacts and focus on certain frequency bands under study. To compare model solutions to experimental data, we apply a Butterworth-bandpass filter of 6th order with lower and higher cut-off frequency 1.25 and 50 Hz, respectively.

RESULTS
This section presents the model topology for certain parameter sets and its dynamics for weak and strong noise, where the noise processes in the networks may be distributed uniformly (iid) or non-uniformly as a superposition of noise processes with different mean values (non-iid). The final application to brain activity indicates the possible importance of the stochastic analysis presented.

Additive iid-Noise
At first, we clarify the topology of Equation (17) providing first insight into the dynamics of the network dynamics. The subsequent numerical simulation of the full model (1) reveals a noise-induced phase transition from a non-coherent state to a phase-coherent state. The complementary analytical description explains this transition.

Equilibria and Linear Stability
The model system (1) exhibits diverse dynamics dependent of the matrices and the external input. From (17), we observe that the means of the network F 0 , M 0 , the maximum firing rate H 0 , and the statistical properties of the external input fully determines the dynamics of the system. To achieve oscillatory activity, we have chosen parameters F 0 , M 0 , I 2 in such a way that bi-stable states may exist for some values of I 1 , D 1 , and D 2 .
Since the additive noise processes {ξ (1) n } are normally distributed, the probability density of υ n at each node n in Equation (19) converges over time to the stationary normal distribution p(υ n ) ∼ N (ξ n , D 1 ). Forξ n = 0 all stationary fluctuations υ n are iid and they obey the same probability density function leading to p V (υ) ∼ N (0, D 1 ) in Equation (14). This yields the transfer function (14)   (24) demonstrating very good accordance. In addition, one observes that larger additive noise level flattens the transfer function and hence affects the nonlinear interactions in the system [37]. The networks V and W excite and inhibit themselves (E − E and I − I interactions) described by the connectivity matrix F with strength F 0 . Their cross-interactions E − I and I − E are given by the matrix M with strength M 0 . In neuroscience, the balance between excitation and inhibition is supposed to play an important role in neural information processing [30,38,39]. In spatially extended biological systems, the balance between excitation and inhibition determines the type of emerging spatiotemporal patterns [40]. To obtain a first insight into the possible dynamics of the system (1), we define the ratio of excitation and inhibition γ = F 0 /M 0 .
Then the equilibrium (V 0 ,W 0 ) of Equation (17) with e 1,2 = 0 is gained by dV/dt = 0, dW/dt = 0. Utilizing Equation (24) and the equivalent formulation for transfer function S 2 and linearization of Equation (17) about the equilibrium the ansatz (V(t),W(t)) = ze λt with the eigenvector z allows to compute the Lyapunov exponents λ 1,2 ∈ C. For complex-valued Lyapunov exponents, i.e., oscillation about the equilibrium, the imaginary part reads where s 1,2 = dS 1,2 (x)/dx are the non-linear gains evaluated at the corresponding equilibrium. Then the eigenfrequency of the oscillation is f = Im(λ)/2π. Figure 2 shows the equilibrium (V 0 ,W 0 ) as a function of the excitation-inhibition ratio γ . The linear stability of the equilibrium is read off the maximum real part of the Lyapunov exponent λ max = max(Re(λ 1 ), Re(λ 2 )). Increasing γ , the system exhibits a single stable focus at low values followed by bistability with top stable node, bottom stable focus, and a center saddle node. For larger γ , the lower stable focus loses stability and transitions to an unstable focus. This transition yields two unstable equilibria and a single top stable node. Further increase of γ causes center and bottom equilibrium to vanish while retaining the top stable node that does not change its stability for larger γ . The frequency of the bottom focus decreases monotonously with increasing γ .
Decreasing the cross-network interaction strength M 0 retains the bi-stability (Figure 3) and M 0 → 0 but makes the lower focus and center node vanish yielding the single equilibrium V 0 → I 1 ,W 0 → I 2 (horizontal line in Figure 3). For large values V 0 ,W 0 the top stable node converges tō see dotted line in Figure 3. Then M 0 → ∞ leads to a single equilibrium at γ = 1/H 0 (vertical dotted line in Figure 3). For M 0 → ∞ but small valuesV 0 ,W 0 , we obtain γ → 1 and S 1 (V 0 ) → S 2 (W 0 ), i.e.,V 0 ≈W 0 ≈ 0. Both latter cases for M 0 → ∞ are illustrated in Figure 3: the larger the cross-interaction strength M 0 the steeper the functionV 0 (γ ) for largeV 0 ,W 0 and the more the right-hand side bistability border stretches to the right. This analysis shows that, for the chosen constant external input and the given noise levels, the bistability between top stable node and bottom oscillatory node is present for γ ≤ 1, i.e., F 0 ≤ M 0 . We point out that different external input may yield bistability for γ > 1 and refer to future work.  In contrast, strong noise (D 1 = 0.8) yields a new lower equilibrium at V 0 n < 0, n = 1, . . . , N (cf. Figures 4A,B). New additive noise-induced states have been studied already close to stability thresholds [9,11,42] for small noise intensities, however not for large additive noise intensities far from stability thresholds that we address here. The new lower state exhibits an augmented spatio-temporal coherence that is observed both as a prominent spectral peak at f ≈ 0.3 Hz in the power spectrum of the spatial averageV(t) and the prominent temporally averaged spatial coherence gPLV (cf. Figure 4C). Similar results are gained for networks with medium mean degree cN but not for low mean degrees (cf. the Appendix).

Numerical Results and Analytical Explanation
Additive noise triggers multistability in the network. Indeed, the analysis of equilibria of Equation (17) with e 1,2 = 0 reveals three states for low noise levels D 1 and a single state for large D 1 (cf. Figure 5). Linear stability analysis of the equilibria (cf. section 3.1.1), further reveals bistability with a stable node in the top branch, an unstable center node and a stable focus at the bottom equilibrium.
The eigenfrequency of the stable focus f diminishes with increasing noise level since the transfer function flattens and hence decreases s 1 (cf. Equation 25). The eigenfrequency approaches f ∼ 0.3Hz for D 1 = 0.8 (cf. Figure 5) and agrees with the oscillation frequency of the spatial mean, shown in Figure 4. Summarizing, this analysis reveals that the transition from a top activity state to a bottom activity state with increasing noise level observed in Figure 4 is a phase transition from a stable node to a stable focus. Figure 4B shows stochastically evolving spatial averagesV. The stochastic force originates from the error of the transfer function ε seen in Figure 1 (cf. Equation 24), and the finite-size  correction terms e 1,2 . Since ε and e 1,2 are random and vary over time, their sum represents a stochastic force in Equation (17).

Extension to Additive Non-iid Noise
We have seen above that the iid additive noise force shapes the non-linear interactions. Now we relax the condition of a vanishing noise average at each node and consider non-iid noise, i.e., a non-Gaussian probability density function p V (υ). Equation (21) allows to describe heterogeneous non-iid noise, i.e., Gaussian processes with diverse meansξ 1 n and variances D (1) n . To demonstrate the effect of non-iid additive noise on our system, we consider M = 2 classes of Gaussian processes with identical variance andξ 1 = −ξ 2 = µ, p 1,2 = 0.5 and distribute the additive noise processes of both classes uniformly over the network. Here µ is the noise mean shift and quantifies the degree of noise heterogeneity. Figure 6A presents the field activity V(t) where µ increases linearly in time. We observe a transition from an upper state to a lower state at µ ≈ 0.68. Figure 5, this phase transition involves a jump from a non-oscillatory non-coherent state to an oscillatory phasecoherent state (cf. Figures 6B,C). To understand this transition, we consider again the averagesV,W and their temporal evolution according to Equation (17). Figures 6D,E show the corresponding equilibria and their linear stability subject to the noise mean shift. For weak noise heterogeneity, i.e., small µ, the system is bistable. Large heterogeneity transforms it to a monostable system with an oscillating state via a saddle-node bifurcation. The eigenfrequency of the oscillation state of f = 0.5Hz is close to the maximum power and phase coherence peak frequency of f ≈ 0.25 Hz in Figures 6B,C. This difference originates from the approximation error of the transfer functions S (cf. Figure 1), resulting to an effective increase in additive noise. For instance, increasing the additive noise variance to D 1 + 0.02 · H 0 F 0 (cf. Figure 1) decreases the eigen frequency at µ = 0.8 in Figure 6 to f = 0.3Hz, i.e., much closer to the observed value. Figure 6F (top panel) compares the analytical (Equation 8) and numerical solution of the probability density p V (υ) for homogeneous and heterogeneous additive noise and confirms good agreement between the analysis and numerics. The corresponding transfer functions show good agreement as well (cf. Figure 6F, bottom panel).

Link to Coherence Resonance
Coherence resonance is present when applied noise induces regular rhythms in a system, although the noise-free system does not [43][44][45]. We have investigated whether the noise-induced transition from non-rhythmic non-coherent system dynamics to rhythmic coherent dynamics shows this phenomenon over a larger range of noise intensities as well. Figure 7 (left) shows the power spectrum ofV(t) subjected to additive iid-noise (cf. Figure 4). We observe no rhythmic activity for small noise  Figure 4; results for W n (t) are equivalent. We have applied the Euler-Maryuama integration technique [41] with time step t = 0.1. intensity (the system evolves on the top branch seen in Figure 4), whereas coherent network activity is present for large noise D 1 ≥ 0.5. In comparison to the analysis shown in Figures 4-7 (left) shows a spectral peak that moves to lower frequencies for larger noise intensities while diminishing in power. Finally, for very large noise intensities D 1 > 10.5, the spectral peak vanishes. This means that coherent network activity is present for intermediate noise levels and the system does not exhibit coherent activity for very small and very large noise intensities. Qualitatively, this resembles coherence resonance.
To evaluate our analytical description, Figure 7 (right) shows the linear power spectrum about the corresponding equilibrium and we observe a qualitatively similar power spectrum subject to the noise intensity, i.e., no phase coherence rhythmic state for very small and very large noise intensities but a spectral peak for intermediate levels. We point out that the analytical power spectrum reproduces quantitatively well the spectral peak frequency at low noise intensities, while the spectral peak already vanishes for noise intensities smaller than the ones observed in the numerical study.
To complete the study, Figure 8 compares the numerical and analytical power spectra similar to Figure 7, but now varying the degree of noise heterogeneity µ. We observe a coherence resonance-like behavior with a maximum coherence at intermediate levels of heterogeneity and vanishing coherence at very low and very high degrees of heterogeneity. Both numerical and analytical spectra are similar qualitatively, but differ in the frequency of maximum coherence. This results from the linear approximation of the analytical power spectrum and the finitesize effect observed and discussed already in the context of

Experimental Results
Finally, to illustrate the importance of the demonstrated results, we consider experimental human brain data. An electroencephalogram (EEG) was recorded during isoflurane induced cardiac anesthesia from a polar FP1 to mastoid derivation at 125 Hz sampling rate using a Narcotrend (R) monitor (MT MonitorTechnik, Hannover, Germany). The EEG was recorded as part of the EPOCAS study (National Clinical Trials number 02976584). A bandpass filter from 0.5 to 45 Hz was applied and cardiac artifacts minimized by subtracting an artificial artifact reference signal from the EEG.
In Figure 9A the upper panel shows an episode with burst suppression, i.e., an intermittent suppression of high-amplitude to small-amplitude activity. Here, burst EEG shows a broad frequency band and the suppressed EEG exhibits oscillations in a narrow frequency band at about 8 Hz (lower panel). This burst suppression was caused experimentally by a rapid increase in isoflurane concentration from 0.8 to 1.8% (not   2) additive noise levels, respectively; circles mark the high and low input with corresponding equilibrium (left). The normalized power spectral density of the burst (blue) and suppression (red) phase shows broad-band and phase-coherent states, respectively (right). In the model simulations, other parameters are taken from Figure 4 and time is re-scaled by t → 0.04t.
shown). This jump from a broad frequency range activity to a coherent oscillation in a narrow frequency band resembles the additive noise effect described in the paragraphs above. The electromagnetic source of EEG lies in the cortex and the reticular activation system (RAS) in the brain sets the excitation level of the cortex [46]. We assume that model (1) describes cortical activity and the additive input is supposed to originate directly or indirectly from the RAS. Under anesthesia, RAS activity varies [47] and hence burst suppression may result from reduced RAS input induced by increased isoflurane concentration. Moreover, synaptic inhibition in local cortical structures may also decrease the intrinsic noise level [36,48]. Hence D 1 reduces for larger isoflurane concentration. Figure 9B (upper panel) shows corresponding simulation results resembling well the experimental EEG. Band-pass filtered simulated activity exhibits a frequency at about 7 Hz during the burst suppression phase and a broad-band activity before and after (cf. Figure 9B, bottom right panel). The underlying model considers global Gaussian spectral-white input I 1 0 (t) = e(I 0 + η(t)), η = 0, η(t)η(τ ) = Aδ(t − τ ) and A is the global noise intensity. For large input the system evolves with a high-amplitude broad band activity about the upper stationary solution and jumps to a small-amplitude narrow band activity about the lower equilibrium ( Figure 9B, lower left panel) for small additive input. We observe that the system is bistable for the large input and monostable for low input. This finding indicates that suddenly decreased and increased input directly or indirectly from the RAS may induce burst suppression in the cortex.

DISCUSSION
In the present work, we have shown that both additive iid noise and non-iid noise, either through the form of fluctuations or heterogeneity, shape the dynamics of our two population model. This occurs because of an effective transformation of the nonlinear interaction function: from a step function for vanishing noise over a symmetric sigmoid for non-vanishing homogeneous noise to a double-sigmoid function for heterogeneous noise. We present a fully analytical description of this system behavior in the context of oscillatory activity, and demonstrate that the system exhibits coherence resonance-like behavior. This enabled us to connect the dynamics to experimental observations, specifically the occurrence of burst suppression in human electroencephalographic data during general anesthesia.
The stochastic analysis implies several major assumptions: • Homogeneity in Equation (2), • Independent and identically distributed stationary stochastic processes in each noise class G m in Equation (4) reflecting ergodicity in each noise class and • Statistical independence of connectivity matrix elements and stochastic process in Equation (12).
These assumptions are major conditions for the mean-field approach considered in the present work and they rule out the application to other network models, such as sparsely connected networks. The network model in the present work is deliberately simple in order to allow the identification of major elements of the novel stochastic analysis. The model assumption of excitationinhibition interaction is motivated by the well-known fact, that such interactions allows oscillatory activity. For instance, the balance between excitation and inhibition is an important mechanism in oscillatory biological neural networks [30,38,39]. Typical cortical neural networks have four times more excitatory than inhibitory synapses. Conversely, the networks V and W have identical number of nodes. This number relation has been chosen for simplicity reasons since the major stochastic effect is independent of the node number. Moreover, balanced networks exhibit an identical level of excitation and inhibition, i.e., γ = 1 in Figures 2, 3. Conversely, parameters in the model, e.g., the constant external inputs and the noise level, have been chosen for simplicity in order to demonstrate and visualize well the stochastic transitions by iid and non-iid noise. These parameters imply γ < 1, i.e., smaller excitation than inhibition, whereas different external input may allow γ > 1 as well. Since a corresponding detailed parameter study would well exceed the major aim of the present work, we refer the reader to forthcoming work.
We point out that the gained results are expected in other network topologies as well, just retaining the homogeneity condition (2) and the statistical independence of connectivity and external noise (12). For instance, previous studies have proposed different network topologies to model the brain that fulfill these conditions [3,49]. However, these conditions may not hold for other models, such as sparsely connected networks [23], clustered networks [38], or scale-free networks [50]. Moreover, we point out that the present model describes the microscopic dynamics of two single interacting populations similar to the Wilson-Cowan model [25] and the homogeneity assumption is reasonable. This means the resulting model describes the dynamics on a larger, mesoscopic, scale. Conversely, previous studies using connectivity data from the brain connectome [3,26] target the global brain dynamics, where each node may obey a Wilson-Cowan population dynamics. In such studies, the underlying connectivity matrices are macroscopic and non-homogeneous.
Our work adds to previous work [13,36,49] highlighting the importance of noise in shaping mesoscopic-scale brain activity. The considered network model may represent a single mesoscopic node in a larger macroscopic network model of the brain. In such macroscopic networks, the noise considered in our model may originate from other mesoscopic populations. Since such distant populations may oscillate and fluctuate randomly, future work will investigate a superposition of oscillatory and noisy external input.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
This study was approved by the Canton of Bern Ethics Committee. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
AH and JL conceived the study and performed all calculations. DH and HK have provided the data and expertise in anesthesia and EEG. All authors have written the manuscript.