Jérémie Lefebvre

Axel Hutt
email: axel.hutt@inria.fr

Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Brain networks are known to generate characteristic, stereotyped activity patterns reflective of self-organization properties. The features of these responses are encoded in space (e.g. as spatial patterns [START_REF] Gardner | Toroidal topology of population activity in grid cells[END_REF] ), in time (e.g. coherent oscillations in given frequency bands [START_REF] González | Communication through coherence by means of crossfrequency coupling[END_REF] ), or a combination of both (e.g. spatio-temporal dynamics such as travelling waves [START_REF] Prechtl | Visual stimuli induce waves of electrical activity in turtle cortex[END_REF] ). Frequently observed experimentally on a mesoscopic scale of few millimeters and frequencies in the range 1-80Hz, such patterns indicates a high degree of self-organisation [START_REF] Haken | Brain Dynamics[END_REF][START_REF] Hutt | Brain connectivity reduction reflects disturbed self-organisation of the brain: Neural disorders and general anaesthesia[END_REF][START_REF] Hutt | Additive noise tunes the self-organization in complex systems[END_REF] . This self-organisation is observed across a wide variety of complex systems [START_REF]Synergetics[END_REF][START_REF] Kaschube | Self-organization and the selection of pinwheel density in visual cortical development[END_REF][START_REF] Kelso | Dynamic Patterns: The Self-Organization of Brain and Behavior[END_REF][START_REF] Nicolis | Self-Organization in Non-Equilibrium Systems: From Dissipative Structures to Order Through Fluctuations[END_REF] , but challenging to characterize mathematically. The specific challenge in neural networks is the emergence of coherent, ordered patterns in spite of the seemingly random, disordered links/connections observed at the microscopic scale [START_REF] Hansel | Modeling feature selectivity in local cortical circuits[END_REF][START_REF] Kaiser | Hierarchy and dynamics of neural networks[END_REF] , spanning scales of tens of micrometers and time scales of few milliseconds. Each neural element, let it be an individual neuron or local population, represents an excitable system with its own dynamics and time scale. Moreover, such elements may further adapt their properties to previous stimuli and display a form of memory. Neuron properties are further known to exhibit a large diversity [START_REF] Rich | Resilience through diversity: Loss of neuronal heterogeneity in epileptogenic human tissue impairs network resilience to sudden changes in synchrony[END_REF] , which determines their stability and resilience to perturbations [START_REF] Hutt | Intrinsic neural diversity quenches the dynamic volatility of neural networks[END_REF] .

In sum, such heterogeneity and seemingly irregular random connectivity structure causes incoming activity to single neurons as part of a network to appear being irregular and random. This holds true , especially in non-sensory brain areas, where neural activity is dominated by recurrent, non-linear interactions and correlates poorly with stimuli [START_REF] Gerstner | Spiking Neuron Models[END_REF] . In such a regime, the relationship linking stimuli and neural responses is blurred and notoriously hard to characterize. Fluctuations in neural activity are hence commonly taken into account through endogenous additive random perturbations. This endogenous noise also takes into account multiple microscopic neurophysiological processes, such as membrane potential fluctuations caused by ion channel gating and/or spontaneous synaptic receptor processes. Such endogenous fluctuations have been shown to have a major impact on neural networks' stability and oscillations, e.g. in the visual cortex [START_REF] Hutt | Arousal fluctuations govern oscillatory transitions between dominant γ and α occipital activity during eyes open/closed conditions[END_REF] , in the corticothalamic feedback system [START_REF] Hutt | Suppression of underlying neuronal fluctuations mediates EEG slowing during general anaesthesia[END_REF] or under various forms of neurostimulation [START_REF] Lefebvre | Stochastic resonance mediates the state-dependent effect of periodic stimulation on cortical alpha oscillations[END_REF] .

Additive input noise may further be used to account for the effect of stimuli and exogeneous inputs. For instance, highfrequency afferent neural signaling or neuromodulatory fluctuations may modulate the noise statistics in a time-dependent way, further influencing neural activity patterns and their properties, such as stability and frequency [START_REF] Solanka | Noise promotes independent control of gamma oscillations and grid firing within recurrent attractor networks[END_REF] . This raises the question whether additive noise -and its impact on neural networks -may be sufficient to explain and reconcile experimental observations. Event-related desynchronization (ERD) and event-related synchronization (ERS) are two important phenomena reliably observed in EEG signals during the processing of sensory, cognitive, or motor events [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF] . ERD refers to a decrease in the power of specific frequency bands in the EEG signal, after the presentation of a task-related event or stimulus. In the α-frequency band (8-12Hz), ERD reflects a state of cortical activation or enhanced neuronal processing associated with the cognitive or motor task [START_REF] Mierau | State-dependent alpha peak frequency shifts: Experimental evidence, potential mechanisms and functional implications[END_REF] . Conversely, ERS in the α-band refers to an increase in the power within the same frequency bands, indicating a state of cortical inhibition or decreased neuronal processing associated with the task [START_REF] Klimesch | EEG alpha oscillations: the inhibition-timing hypothesis[END_REF] . ERD/ERS also occurs in other frequency bands. Figure 1 shows experimental source activity in sensorimotor areas observed during walking. One can observe temporal enhancement and reduction of spectral power within the γ-band (25 -50Hz), signatures of ERD/ERS resulting from fluctuating arousal and associated recruitment of neural populations. In general, both ERD and ERS provide valuable insights into the neural mechanisms underlying different cognitive processes and motor actions, observed during cognitive tasks as event-related potentials [START_REF] Iijima | Eventrelated synchronization and desynchronization of high-frequency electroencephalographic activity during a visual Go/No-Go paradigm[END_REF] or during neurostimulation [START_REF] Wilson | tDCS modulates visual gamma oscillations and basal alpha activity in occipital cortices: Evidence from MEG[END_REF] . We hypothesise that both phenomena may be understood mathematically by investigating the response of neural networks to various degrees of additive noise. This is in line with previous studies indicating that endogenous noise may tune brain rhythms [START_REF] Hutt | Suppression of underlying neuronal fluctuations mediates EEG slowing during general anaesthesia[END_REF][START_REF] Solanka | Noise promotes independent control of gamma oscillations and grid firing within recurrent attractor networks[END_REF] and thus resembles coherence resonance [START_REF] Hutt | Coherence resonance in random erdos-renyi neural networks : mean-field theory[END_REF][START_REF] Pikovsky | Coherence resonance in a noise-driven excitable system[END_REF] . In other words, power modulations as in Fig. 1 are supposed to originate from in underlying modulations of endogenous noise levels.

Previous studies indicate that α-rhythms originate from the feedback circuit involving cortical and thalamic areas [START_REF] Halgren | The generation and propagation of the human alpha rhythm[END_REF][START_REF] Hashemi | How the cortico-thalamic feedback affects the EEG power spectrum over frontal and occipital regions during propofol-induced anaesthetic sedation[END_REF][START_REF] Robinson | Modal analysis of corticothalamic dynamics, electroencephalographic spectra and evoked potentials[END_REF][START_REF] Vijayan | Thalamic model of awake alpha oscillations and implications for stimulus processing[END_REF] . Conversely, there is strong evidence that γ-rhythms originate from local circuits [START_REF] Ray | Different origins of gamma rhythm and highgamma activity in macaque visual cortex[END_REF][START_REF] Whittington | Multiple origins of the cortical gamma rhythm[END_REF] . Previous studies showed that a rather simple neural mass model permits to describe γ-rhythms, whereas a much more complex neural mass model is required to model α-activity [START_REF] Hutt | Arousal fluctuations govern oscillatory transitions between dominant γ and α occipital activity during eyes open/closed conditions[END_REF][START_REF] Riedinger | Mathematical model insights into eeg origin under transcranial direct current stimulation (tdcs) in the context of psychosis[END_REF] . For simplicity, the present work employs the simple model to demonstrate how noise may tune γ-rhythms. Here we propose that spectral power modulations (i.e., ERD/ERS) in the γ-frequency range 30 -60Hz (cf. Fig. 1) occur naturally through fluctuations in the amplitude of endogenous neural noise. Such noise-driven amplitude modulation reflects a synchronisation (power enhancement; ERS) or desynchronisation (power reduction; ERD) of neural populations [START_REF] Haig | Synchronous cortical gamma-band activity in task-relevant cognition[END_REF][START_REF] Iijima | Eventrelated synchronization and desynchronization of high-frequency electroencephalographic activity during a visual Go/No-Go paradigm[END_REF][START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF][START_REF] Shibata | Event-related dynamics of the gamma-band oscillation in the human brain: information processing during a Go/No-Go hand movement task[END_REF] . To reduce the complexity of corresponding random network models, meanfield approaches provide important insights. We thus implemented coupled Erdös-Renyi networks of finite number of excitatory and inhibitory nodes and derive a corresponding mean-field model, which describes the temporal dynamics of the network average. Previous theoretical studies on similar networks [START_REF] Hutt | Additive noise-induced system evolution (ANISE)[END_REF][START_REF] Hutt | Arousal fluctuations govern oscillatory transitions between dominant γ and α occipital activity during eyes open/closed conditions[END_REF][START_REF] Hutt | Phase coherence induced by additive gaussian and non-gaussian noise in excitable networks with application to burst suppression-like brain signals[END_REF][START_REF] Hutt | Coherence resonance in random erdos-renyi neural networks : mean-field theory[END_REF][START_REF] Rich | Neurostimulation stabilizes spiking neural networks by disrupting seizure-like oscillatory transitions[END_REF] have considered primarily the case of infinite size, whereas finite-size networks may introduce additional random effects [START_REF] Hutt | Phase coherence induced by additive gaussian and non-gaussian noise in excitable networks with application to burst suppression-like brain signals[END_REF][START_REF] Hutt | Coherence resonance in random erdos-renyi neural networks : mean-field theory[END_REF] which may play an important role in ERD/ERS phenomena. These studies have shown that the mean and variance of additive input noise enters the meanfield equation as control parameter, tuning the stability and spectral properties of the system. If the mean and variance of exogenous noise input changes over time, the mean-field is driven parametrically and may exhibit nonlinear effects, such as stochastic resonance [START_REF] Lefebvre | Stochastic resonance mediates the state-dependent effect of periodic stimulation on cortical alpha oscillations[END_REF] or coherence resonance [START_REF] Baspinar | Coherence resonance in neuronal populations: mean-field versus network model[END_REF][START_REF] Hutt | Coherence resonance in random erdos-renyi neural networks : mean-field theory[END_REF] . The present work will explore in more detail the finite-size impact in random networks. Collectively, our results demonstrate that both synchronisation and desynchronisation in the γ-frequency range can be explained by modulation in noise variance.

The subsequent section introduces the microscopic network under study and shows the derivation of a corresponding stochastic mean-field model of finite size. In the following section, simulation studies demonstrate different types of mean-field dynamics dependent on how noise input is modulated.

II. THE NETWORK MODEL

The present work considers a large random network of number of nodes N. Each node represents an ensemble of microscopic neurons, a so-called micro-column. The network of micro-columns comprises two interacting sub-networks of excitatory and inhibitory nodes or cells. Mathematically, the excitatory and inhibitory nodes are collected in the vectors V ∈ R N and W ∈ R N , respectively, which represent mean dendritic currents of the corresponding micro-column (also known as minicolumn [START_REF] Bennett | An attempt at a unified theory of the neocortical microcircuit in sensory cortex[END_REF][START_REF] Buxhoeveden | The minicolumn hypothesis in neuroscience[END_REF][START_REF] Mountcastle | Modality and topographic properties of single neurons of cat's somatic sensory cortex[END_REF] ). The associated dynamics of these micro-columns obey

τ e dV = (-V + FS 1 [V] -MS 2 [W] + eI e )dt + dY e (t) τ i dW = (-W + MS 1 [V] -FS 2 [W] + eI i )dt + dY i (t) (1)
with the N-dimensional identity vector e = (1, . . . , 1) , the intra-subnetwork random connectivity matrix F , the intersubnetwork random connectivity matrix M and the synaptic time scales τ e and τ i . The connectivity matrices F and M are N × N-dimensional random matrices. Both matrices share the structure of a random matrix K ∈ R N×N , whose elements {(K) i j } are Bernoulli-distributed with probability c and are equal to (K) i j = XK o /cN, K o > 0, where X is a Bernoulli random variable . Then the mean connectivities are

µ F = F 0 /N, µ M = M 0 /N with respective variance σ 2 F = F 2 o (1 -c)/cN 2 , σ 2 M = M 2 o (1 -c)/cN 2 . The external stochastic input is a N-dimensional Wiener pro- cess Y e,i ∈ R N , for which (dY e,i ) n (t) = 0 (dY e,i ) n (t)(dY e,i ) m (t ) = 2D e,i δ nm δ (t -t )dtdt
for n, m = 1, . . . , N and the noise variance D e,i . The symbol

• denotes the ensemble average, δ nm is the Kronecker symbol and δ (•) denotes the δ -distribution.

Nodes in the network exhibit a nonlinear transfer function S 1 (•) and S 2 (•) for excitatory and inhibitory activity, respectively. Such a transfer function at a node represents the population firing rate in the respective micro-column. We assume a

McCullogh-Pitts neuron model [START_REF] Abeles | Corticonics[END_REF][START_REF] Hutt | Study of GABAergic extra-synaptic tonic inhibition in single neurons and neural populations by traversing neural scales: application to propofol-induced anaesthesia[END_REF] with S 1 (x) = S 0 Θ(x), S 0 > 1 and S 2 (x) = Θ(x) with the Heaviside function

Θ(x) = 0 ∀x < 0, Θ(x) = 1 ∀x ≥ 0. The function vector S 1,2 [x] ∈ R N is de- fined by (S 1,2 [x]) n = S 1,2 (x n
). This neuron model implies identical threshold for all micro-columns.

Mean-field theory

The key idea of the derivation is a projection of the network activity on the matrices' eigenmodes, yielding a separation between mean-field-and stochastic modes. The Proposition 1 (see Appendix), valid for both matrices F and M, will prove helpful. Assuming the basis composition

V(t) = N ∑ n=1 x n (t)φ e n (2) W(t) = N ∑ n=1 y n (t)φ i n with the vector basis {φ e n }, {φ i n } and φ e n , φ i n ∈ C N and a bi- orthonormal basis {ψ e k }, {ψ i k } and ψ e k , ψ i k ∈ C N with ψ e, † k φ e l = δ kl , ψ i, † k φ i l = δ kl , ψ x, † k φ y l = 0 ∀ x, y = e, i , Proposition 1 stipulates φ e 1 = φ i 1 = e , ψ e 1 = ψ i 1 = e/N
. Moreover, we define the eigenvalue problems

ψ e,i, † k F = λ F k ψ e,i, † k , ψ e,i, † k M = λ M k ψ e,i, † k (3)
for k = 1, . . . , N.

Then, projections on the vector basis and using Eqs. (3) together with Proposition 1 yields

τ e dx 1 = (-x 1 + λ F 1 N e S 1 [V] - λ M 1 N e S 2 [W] + I e )dt + dη e (t)
τ i dy 1 = (-y 1 + λ M 1 N e S 1 [V] - λ F 1 N e S 2 [W] + I i )dt + dη i (t) (4)
with dη e,i (t) = e dY e,i (t)/N and the remaining modes obey

τ e dx k (t) = (-x k (t) + λ F k ψ e † k S 1 [V] -λ M k ψ e † k S 2 [W] + ψ e † k eI e )dt + dη e k (t) τ i dy k (t) = (-y k (t) + λ M k ψ i † k S 1 [V] -λ F k ψ i † k S 2 [W] + ψ i † k eI i )dt + dη i k (t) (5)
for k > 1 with dη e,i k (t) = ψ e,i, † k dY e,i (t). The mean-field modes obey Eqs. ( 4) and determine the stochastic modes in Eqs. (5), which in turn influences the mean-field modes' evolution.

The Proposition 2 (see Appendix) details how both groups of modes are interrelated, how the implicit closure problem in mean-field theory is solved in the current context, and how the stochastic modes tune the nonlinear interactions in the network. For the case of the excitatory population, one may rewrite V(t) in Eq. ( 2) as

V = (x 1 -Z(t) +U F S 1 [x 1 ] -U M S 2 [y 1 ])e + w e (t)
with

w e (t) = 1 τ t -∞
e -(t-t )/τ dY e (t )

Z(t) = 1 τ t -∞
e -(t-t )/τ dR e (t ) .

and the finite size random corrections

U F ∼ N (0, σ 2 F /N), U M ∼ N (0, σ 2 M /N).
The finite size-fluctuations R e (t) are Wiener processes with dR e ∼ N (0, D e /N)dt, the finite size-fluctuations Z(t) obeys an Ornstein-Uhlenbeck process or in physical terms red-coloured noise.

Moreover, assuming the function vector

S[u] = (S 1 (u 1 ), . . . , S 1 (u N )) ∈ R N like in Eq. (1) and u = xe + w ∈ R N with x ∈ R and random vector w ∈ R N , w i ∼ N (0, D e /τ e ), Proposition 2 yields e N S[u] = S(x) + S 2 0 S(x) -S2 (x) N χ(t)
with random variable χ(t) ∼ N (0, 1) and

S(x) = ∞ -∞ S 1 (x + w)p e (w)dw = S 0 2 1 + erf x 2(D e /τ e ) .
The dynamics of W(t) can also be written similarly and the derivation resembles closely to the one of V(t).

The application of Proposition 2 is heavily based on the homogeneity assumption of the network: the network degree at each node is very similar across nodes, and identical in infinite networks, i.e. the degree distribution is unimodal with a narrow peak. This results into a homogeneous eigenvector in the matrices' edge spectrum i.e., φ e 1 ≈ e, and consequently the impact of the transfer function is identical at each node. Hence, the nonlinear interaction and the network decouple. If the network topology does not exhibit such a narrow degree distribution (e.g. in scale-free networks or small-world networks), the node connectivity is not identical in all nodes. Hence, the eigenvector becomes φ e 1 = e and Eq. ( 4) does not hold anymore. This implies that the impact of the nonlinear transfer function is different at each node and hence nonlinear interactions and network do not decouple. The results of Proposition 2 allow us to recast Eqs. (4) as

τ e dx 1 = -x 1 + λ F 1 G e (x 1 -Z(t) + T 1 (x 1 , y 1 ),t) -λ M 1 G i (y 1 -W (t) + T 2 (x 1 , y 1 ),t) + I e dt + dR e (t) τ i dy 1 = -y 1 + λ M 1 G e (x 1 -Z(t) + T 1 (x 1 , y 1 ),t) -λ F 1 G i (y 1 -W (t) + T 2 (x 1 , y 1 ),t) + I i dt + dR i (t) (6)
with

G e (x 1 ,t) = Ḡe (x 1 ) + ∆G e (x 1 ) • X e (t) G i (y 1 ,t) = Ḡi (y 1 ) + ∆G i (y 1 ) • X i (t) , (7)
for temporal fluctuations X e,i (t) ∼ N (0, 1) and where

Ḡe (x 1 ) = ∞ -∞ S 1 (x 1 + w) p e (w)dw Ḡi (y 1 ) = ∞ -∞ S 2 (y 1 + w) p i (w)dw ∆G e (x 1 ) = S 2 0 Ḡe (x 1 ) -Ḡ2 e (x 1 ) N ∆G i (y 1 ) = Ḡi (y 1 ) -Ḡ2 i (y 1 ) N .
Combining previous results, the network exhibits finite-size fluctuations dR e = ξ e dt , ξ e ∼ N (0, D e /N)

dR i = ξ i dt , ξ i ∼ N (0, D i /N) Z(t) = 1 τ t -∞
e -(t-t )/τ dR e (t )

W (t) = 1 τ t -∞
e -(t-t )/τ dR i (t ) ,

where random state-dependent corrections of the nonlinear interactions occur through the terms ∆G e,i , which now possess random state-dependent thresholds

T 1 (x 1 , y 1 ) = U F G e (x 1 ) -U M G i (y 1 ) T 2 (x 1 , y 1 ) = U M G e (x 1 ) -U F G i (y 1 ) , with U F ∼ N (0, F 2 0 (1 -c)/cN 3 ), U M ∼ N (0, M 2 0 (1 - c)/cN 3 ).
For large networks,

F 2 0 (1 -c)/cN 3 D e,i /N, M 2 0 (1 -c)/cN 3 D e,i
/N and thus the statedependent thresholds are negligible compared to the finite-size fluctuations. If N → ∞, T 1,2 → 0, i.e. these terms are finite-size threshold corrections and vanish for infinite networks.

Since w(t) obeys an Ornstein-Uhlenbeck process, its stationary probability density function may be written as p e,i (w) = 1 2πD e,i /τ e,i e -w 2 /2(D e,i /τ e,i )

and consequently, according to Proposition 2,

Ḡe (x) = S 0 2 1 + erf x 2(D e /τ e ) Ḡi (x) = 1 2 1 + erf x 2(D i /τ i ) . (8)
III. SIMULATIONS

In the following, we consider large but finite-size networks obeying Eqs. ( 6) with T 1,2 ≈ 0. These latter terms result from the random connectivity structure, but are negligible for large networks as shown above. Then

dx 1 = - x 1 τ e + λ F 1 τ e G e (x 1 -Z(t),t) - λ M 1 τ e G i (y 1 -W (t),t) + J e 1 τ e dt + 1
τ e dR e (t)

dy 1 = - y 1 τ i + λ M 1 τ i G e (x 1 -Z(t),t) - λ F 1 τ i G i (y 1 -W (t),t) + J i 1 τ i dt + 1 τ i dR i (t) (9)
Please recall that G e,i are state-dependent finite-size fluctuations and Z(t), W (t) fluctuate in time according to an Ornstein-Uhlenbeck process with network-size dependent noise strength. The fluctuation strength of Z(t) and W (t) are proportional to √ D e and √ D i , respectively, while the nonlinear gain of the transfer functions G e and G i (their slope) decrease nonlinearly with √ D e and √ D i , respectively. Consequently, stronger additive noise enhances the finite-size fluctuations Z(t),W (t) but diminishes the response sensitivity (the nonlinear gain) of the system in accordance to previous studies [START_REF] Hutt | Additive noise tunes the self-organization in complex systems[END_REF][START_REF] Hutt | Dynamic control of synchronous activity in networks of spiking neurons[END_REF] .

A. Synchronisation modulated by noise

Experimental studies [START_REF] Haig | Synchronous cortical gamma-band activity in task-relevant cognition[END_REF][START_REF] Iijima | Eventrelated synchronization and desynchronization of high-frequency electroencephalographic activity during a visual Go/No-Go paradigm[END_REF][START_REF] Shibata | Event-related dynamics of the gamma-band oscillation in the human brain: information processing during a Go/No-Go hand movement task[END_REF] have revealed transient spectral power of γ-rhythms (frequency range [25Hz; 60Hz]) mirroring cognitive task demands. These are unequivocal manifestations of ERD/ERS. Changes in endogeneous noise amplitude can be shown to generate such effects. Indeed, since microscopic desynchronized (resp. synchronized) activity yields low (resp. high) amplitudes mean-field activity, experimentally observed spectral power amplitude reflects the (changing) degree of synchronisation amongst neural networks. Such synchronization, and how it depends on noise amplitude, can be assessed through mean-field approaches.

ERD by increasing endogenous noise : passing the Hopf bifurcation

For N = 100, Fig. 2(A) shows the network activity of V(t) (upper panel) over time for changing noise variances for the excitatory population (i.e.,D e ) and the corresponding timefrequency plot of the numerical mean V (t) = ∑ N i=1 V i (t)/N (lower panel). Throughout, the inhibitory noise variance (i.e. D i ) is kept constant. For low noise variances, we observe a mean activity V (t) in the γ-frequency range fluctuating in amplitude over time. Increasing the noise variance suppresses this γ-activity. Returning to the lower noise variance, γ-activity re-emerges again with fluctuating amplitude. This behaviour can be understood by a closer look at the stationary state x1 resulting from the fixed point condition dx 1 /dt = dy 1 /dt = 0 in the absence of finite-size noise, i.e. N → ∞. Linearizing Eqs. (9) about the stationary states x1 , ȳ1 with x 1 (t) = x1 + u(t), y 1 = ȳ1 + v(t) we obtain the following dynamical system

τ dz/dt = Lz + s(t) (10)
with z = (u, v) and

τ = τ e 0 0 τ i , L = λ F 1 G e -1 -λ M 1 G i λ M 1 G e -λ F 1 G i -1 (11) s = -λ F 1 G e Z(t) + λ M 1 G i W (t) -λ M 1 G e Z(t) + λ F 1 G i W (t) . (12)
Here, G e,i = dG e,i (x)/dx is evaluated at x = x1 , ȳ1 . The stationary states' stability is determined by the maximum eigenvalue real part λ m of L. If the eigenvalue λ m is complexvalued with I m(λ m ) = 2πν, then ν is called eigenfrequency. (B) Instantaneous power in the γ-frequency range. A Butterworth online filter of 4th order was applied to the mean-field activity, the filter output was squared and averaged in a sliding window of 2s duration. We observe an increase in power reflecting enhanced synchronisation in the network. Other parameters were identical to parameters used in Fig. 2. The influence of noise on oscillatory activity can be understood through the non-linear mechanism of hysteresis. For low noise, three stationary states exist: the upper branch is a stable node, the centre branch is a saddle node and the bottom branch is an unstable focus, cf. the corresponding maximum eigenvalue real part λ m shown in the lower panel. We observe a Hopf-bifurcation on the lower stationary state. Increasing the excitatory noise variance D e moves the system over the Hopf bifurcation point and thus stabilises the focus and merges the stable node and the saddle node by a saddlenode bifurcation. Here, the Hop-bifurcation and the saddlenode bifurcation occur on different stationary state branches and at different noise levels. Moreover, the eigenfrequency ν of the oscillation on the lower branch decreases monotonically with increasing noise level. These properties explain the simulation results presented in Figure 2(A): the γ-activity observed at low noise levels represents a nonlinear limit cycle above the Hopf-bifurcation point, while the low-amplitude γ-activity at larger noise levels represents a quasi-cycle about the stable focus below the Hopf-bifurcation point.

ERD by decreasing endogenous noise: phase transition

The topology of the mean-field differential equation system (9) provides a justification to a decrease of γ-activity triggered by noise modulation, cf. Fig. 3. If the noise level first decreases at the lower stationary state branch, the system passes the Hopf bifurcation point yielding increasing oscillation amplitude. Then a nonlinear jump from the lower to the . Finite-size impact on oscillation frequency. The spectral power of the quasi-cycle on the lower branch shown in Fig. 2(B) at D e /τ e = 0.4. The power distribution has been computed applying the Bartlett-Welch method for simulated signals V (t) = ∑ N i=1 V i (t)/N (black) as solution of Eq. ( 1) and x 1 (t) (red) as solution of Eqs. (9). Other parameters were identical to parameters used in Fig. 2. upper stationary state branch yields a non-oscillatory state and the γ-amplitude vanishes. Increasing the noise level again induces a second jump from the upper to the lower branch yielding γ-rhythms again. Hence noise reduction induces γ-desynchronisation.

ERS by decreasing noise level: passing the Hopf bifurcation

Similar to the γ-desynchronisation described above, synchronisation may also occur while modulating the noise level in the other direction. Typically, experimental studies investigating synchronisation in evoked potentials apply a frequency band filter and compute the time-dependent power value in this frequency band [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF] . To better compare our results to such an analysis, Fig. 4 presents the time-dependent power in the γ-frequency band while abruptly decreasing the noise level. We observe a transient increase and decrease in γ-power for a noise level decrease and increase, respectively.

B. Finite-size impact

Finally, we consider the impact of finite-size fluctuations in Eq. ( 9) on the results above. Figure 5 shows the power spectrum on the lower stationary state branch at a fixed noise level at different network sizes N. We observe that small networks exhibit a broader spectrum of rhythms about lower fre- quencies compared to large networks. Hence the network size determines the system's principal frequencies and their distribution.

ΔG(V) ΔG'(V) 0 0,5 -0,5 V -1 0 1 G(V) G'(V) 0 0,5 1,0 V -1 0 1 (A) (B)
To explain this, let us take closer look at the finite-size contributions to the nonlinear transfer function G e,i in Eq. ( 9). The transfer functions G e,i in Eqs. ( 7), ( 8) have the form

G(V ) = Ḡ(V ) + ∆G(V ) • X Ḡ(V ) = 1 2 1 + erf V √ 2σ ∆G(V ) = Ḡ(V ) -Ḡ2 (V ) N (13)
with the random variable X ∼ N (0, 1) and real parameter σ > 0. Since the system's frequency distribution depends on the nonlinear gain G e,i , cf. Eqs. ( 10)-( 12), finite-size effects may translate to their modulation and we find

G (V ) = Ḡ (V )       1 + 1 -2 Ḡ(V ) 2 Ḡ(V ) -Ḡ2 (V )N =∆G (V ) •X       . (14)
by taking the derivative of G(V ) in Eq. ( 13). Figure 6(A) presents the finite-size correction to the transfer function ∆G from Eq. ( 13) and its derivative ∆G from Eq. ( 14). The absolute value of the nonlinear gain finite-size correction is maximum about V = ±0.5, which coincides with the approximate lower stationary state location. This indicates that the finite network size affects primarily the state evolving in the γ-frequency range. Figure 6(B) illustrates the impact of fluctuations on the transfer function and its nonlinear gain and we observe larger fluctuations about V = ±0.5 in accordance to (A).

IV. DISCUSSION

Endogeneous noise drives the neural network under study and affects its stability and oscillation rhythms.

The stochastic network dynamics

We find that the network activity is the sum of a mean-field activity, background Ornstein-Uhlenbeck noise processes and finite-size fluctuations, cf. Eq. (B5). In general, the meanfield dynamics depends on both the background noise process and the finite-size fluctuations. For very large networks, finite-size fluctuations are negligible and the mean-field activity obeys a deterministic dynamics, which still depends implicitly of the endogenous noise variance driving the network and hence the background noise process. An analytical study of this deterministic model reveals the networks stationary states, their linear stability and the systems characteristic frequencies. This mathematical insight permits to derive conditions under which the network exhibits coherent and noncoherent behavior subjected to the endogenous noise level and thus exhibits synchronisation and desynchronisation, respectively (Figs. 234). Conversely, smaller networks induce finitesize fluctuations which both couple into-and add up to the mean-field activity. Our studies show that they impact the systems frequency range (Fig. 5) resulting from finite-size effects on the systems nonlinear transfer functions.

ERS and ERD as noise-driven phenomena

Modulating the endogenous noise level permits to tune the system's oscillatory activity and enhance or diminish the magnitude of γ-activity activity and thus induce ERS and ERD, respectively. ERD may be evoked by excitatory noise enhancement by passing a Hopf-bifurcation or it may result from an excitatory noise level decrease evoking a phase transition from an oscillatory to an non-oscillatory state. Conversely, ERS can be explained by reducing the excitatory noise level while passing a Hopf-bifurcation. Experimental studies on motor behavior [START_REF] Andrew | Event-related coherence as a tool for studying dynamic interaction of brain regions[END_REF][START_REF] Bressler | Episodic multi-regional cortical coherence at multiple frequencies during visual task performance[END_REF] indicate that ERD and ERS reflects the activity between and during motor action, respectively. In terms of our model, this indicates alternating noise levels during motor behavior in accordance with previous studies on the role of endogenous fluctuations [START_REF] Boly | Baseline brain activity fluctuations predict somatosensory perception in humans[END_REF][START_REF] Fox | Intrinsic fluctuations within cortical systems account for intertrial variability in human behavior[END_REF] . Moreover, previous other work has provided evidence that ERD may result from reduced network inhibition [START_REF] Takemi | Event-related desynchronization reflects downregulation of intracortical inhibition in human primary motor cortex[END_REF] . This work suggests to investigate in our model whether ERD may be induced by tuning the maximum of transfer function S 2 in Eq. (1). Since this additional study would exceed the major aim of the current study, we refer the reader to forthcoming work.

Limitations and outlook

The network model assumes simple McCullough-Pitts neurons, i.e. threshold units, whereas biological neurons exhibit a much more complex dynamics. Future work should extend the presented network setup to such biologically more realistic models. However, we expect a similar network behavior since the resulting mean-field models will include similar sigmoidal-shape transfer functions.

Moreover, the oversimplifying assumption of an Erdös-Renyi topology can be relaxed easily and extended to Gaussian random networks. Both network types assume a homogeneous network with a narrow uni-modal degree distribution. Proposition 2 is important in the mean-field derivation and applies only if this homogeneity assumption holds and the nonlinear interaction and network decouple . Extensions to more biologically realistic topologies, such as scale-free [START_REF] Holme | Rare and everywhere: Perspectives on scale-free networks[END_REF][START_REF] Lee | Brain networks maintain a scale-free organization across consciousness, anesthesia, and recovery: evidence for adaptive reconfiguration[END_REF] or smallworld networks [START_REF] Bassett | Small-world brain networks revisited[END_REF] will be more challenging due to the resulting coupling of network and transfer function. For corresponding more detailed discussions, we refer the reader to forthcoming work. Similarly, more realistic network models exhibit a spatial embedding yielding space-time dynamics observed in neural tissue [START_REF] Castelnovo | Sleep spindles and slow waves in schizophrenia and related disorders: main findings, challenges and future perspectives[END_REF][START_REF] Huang | Spiral waves in disinhibited mammalian neocortex[END_REF] . The present work assumes additive noise input, which is uncorrelated in space and time. More realistic input may exhibit correlations, which may influence the network response differently [START_REF] Churchland | Stimulus onset quenches neural variability: a widespread cortical phenomenon[END_REF][START_REF] Zhou | Impact of neuronal heterogeneity on correlated colored noise-induced synchronization[END_REF] . Moreover, the present work does not investigate how finite-size fluctuations interact to shape spectral power. Previous studies indicate an important role of such fluctuations in neural processing [START_REF] Pietras | Mesoscopic description of hippocampal replay and metastability in spiking neural networks with shortterm plasticity[END_REF] . and σ 2 K , σ 2 1 → 0 for N → ∞, it follows that (φ 1 ) i → 1. This last statement is based on the corresponding results for symmetric Bernoulli matrices for very large networks with cN → ∞ 13 only, but is confirmed by numerical verification.

The mean eigenvalue of the bulk spectrum is E[λ k ] = 0 with maximum value [START_REF] May | Stability and complexity in model ecosystems[END_REF][START_REF] Wigner | On the distribution of the roots of certain symmetric matrices[END_REF]

λ k>1,max = 2σ K √ N (B3) = 2K 0 1 -c cN (B4)
and the entries of their eigenvectors are random [START_REF] Ding | Spectral distributions of adjacency and laplacian matrices of random graphs[END_REF][START_REF] Guilonnet | Bernoulli random matrices[END_REF] .

We observe that |λ k>1,max | → 0 for N → ∞.

2. For the orthonormal basis {ψ}, we assume

ψ † k K = λ k ψ † k
and we obtain

Kφ k = λ k φ k ψ † k Kφ k = λ k ψ † k φ k ψ † k Kφ k = λ k λ k ψ † k φ k = λ k λ k = λ k .
Since φ 1 = e, we obtain ψ 1 = e/N by relation (B1).

Proposition 2. The relation of Eqs. (4) and Eqs. ( 5) is described in the first part, after which the nonlinear interaction is detailed in the second part.

1. The mode expansion (2) can be re-written as

V = (x 1 -Z(t) +U F S 1 [x 1 ] -U M S 2 [y 1 ])e + w e (t) (B5)
with

w e (t) = 1 τ t -∞
e -(t-t )/τ dY e (t )

Z(t) = 1 τ t -∞
e -(t-t )/τ dR e (t ) .

and the finite size random corrections U

F ∼ N (0, σ 2 F /N), U M ∼ N (0, σ 2 M /N).
The finite size-fluctuations R e (t) are Wiener processes with dR e ∼ N (0, D e /N)dt, the finite size-fluctuations Z(t) obeys an Ornstein-Uhlenbeck process, or in physical terms, red-coloured noise.

2. Further assume the function vector S[u] = (S 1 (u 1 ), . . . , S 1 (u N )) ∈ R N like in Eq. ( 1) and u = xe + w ∈ R N with x ∈ R and random vector w ∈ R N , w i ∼ N (0, D e /τ e ). Then

e N S[u] = S(x) + S 2 0 S(x) -S2 (x) N χ(t) (B6)
with random variable χ(t) ∼ N (0, 1) and

S(x) = ∞ -∞ S 1 (x + w)p e (w)dw = S 0 2 1 + erf x 2(D e /τ e ) .
Proof.

1. We consider

τ e dx k (t) = (-x k (t) + ψ e † k Q k )dt + dη e k (t) (B7)
for k = 2, . . . , N motivated by Eq. ( 5) with

Q k (t) = λ F k S 1 [V(t)] -λ M k S 2 [W(t)] + eI e .
The solution of (B7) is

x k (t) = t -∞ e -(t-t )/τ e dη e k (t ) + t -∞ e -(t-t )/τ e ψ e † k Q k (t )dt
and Eq. ( 2) reads then

V = x 1 φ e 1 + N ∑ k=2 φ e k t -∞
e -(t-t )/τ e ψ e, † k dY e (t )

+ N ∑ k=2 t -∞ e -(t-t )/τ e ψ e † k Q k (t )dt . ( B8
)
This equation is implicit in V. This is a manifestation of the closure problem, well-known in mean-field theory. The problem is solvable uniquely in specific cases under certain assumptions. Here, we assume that V(t) evolves very slowly over time, i.e. slower than the time scale τ e and thus Q(t) ≈ const. This implies that Q(t) does not fluctuate randomly. Consequently,

Q k = λ F k S 1 [ V] -λ M k S 2 [ W]
+ eI e and V ≈ V, W ≈ W is constant in time. Then the mode projection amplitudes x k (t) obey Ornstein-Uhlenbeck processes

x k (t) = with w i ∼ N (0, D e /τ e ).

Figure 1 .

 1 Figure 1. Synchronisation (ERS) and Desynchronisation (ERD) in human experimental source activity in sensorimotor areas during walking. The time-frequency distribution of relative spectral power change in trial-averaged single-subject data reveals a sequence of alternating ERS (red) and ERD (blue) in the γ-frequency band [30Hz; 50Hz]. Reproduced from Fig. 3 in Seeber et al. 56 licensed under a Creative Commons Attribution (CC BY) license.

Figure 2 .

 2 Figure 2. γ-desynchronisation by increasing excitatory input noise (A) Top panel: network activity V i obeying Eq. (1) with the initial state close to the lowest stationary state. The initial noise variance was chosen to D e = 0.25τ e and changed to D e = 0.8τ e at t = 5s and changed back to D e = 0.25τ e at t = 15s. Bottom panel: Time-frequency plot of instantaneous spectral power (windowed Fourier Transform connvolved with a Hamming window with 2s window width and overlap of 1.8s). (B) Top panel: the stationary state x1 dependent on D e /τ e revealing s stable node (dotted line), saddle node (dashed line) and a stable focus (solid line). Bottom panel: the maximum eigenvalue real part λ m (left hand side) and the eigenfrequency of the focus (right hand side). Other parameters were D i /τ i = 0.2, N = 100, F 0 = 2.17, M 0 = 3.87, τ e = 0.005, τ i = 0.02, S 0 = 1.7 . and the integration time step was 0.5ms utilising the Euler-Maruyama method.

Figure 3 .

 3 Figure 3. γ-desynchronisation by decreasing excitatory input noise (A) Top panel: the time-dependent change of the noise variance D e /τ e in the simulation. Bottom panel: x 1 (t) simulated from Eq. (9) dependent on D e /τ e and thus on time. It shows a jump from the lower stationary state branch to the upper branch and back. (B) The time-dependent change of D e /τ e moves the system's stationary state x1 and induces noise-induced jumps between the lower focus and the upper stable node, cf. (A), bottom panel. The jumps from the bottom to top branch and back occur at different values of D e /τ e reflecting a hysteresis. The period on the top branch reflects a nonoscillatory, i.e. desynchronised, behavior. (C) Time-frequency plot of spectral power of x 1 (t). Other parameters were identical to parameters used in Fig. 2.

Figure 4 .

 4 Figure 4. γ-synchronisation by decreasing excitatory input noise (A) Top panel: the time-dependent relation D e /τ e in the simulation.(B) Instantaneous power in the γ-frequency range. A Butterworth online filter of 4th order was applied to the mean-field activity, the filter output was squared and averaged in a sliding window of 2s duration. We observe an increase in power reflecting enhanced synchronisation in the network. Other parameters were identical to parameters used in Fig.2.

Figure 2 (

 2 Figure 2(B,top panel) shows the stationary state x1 for different excitatory noise levels D e .The influence of noise on oscillatory activity can be understood through the non-linear mechanism of hysteresis. For low noise, three stationary states exist: the upper branch is a stable node, the centre branch is a saddle node and the bottom branch is an unstable focus, cf. the corresponding maximum eigenvalue real part λ m shown in the lower panel. We observe a Hopf-bifurcation on the lower stationary state. Increasing the excitatory noise variance D e moves the system over the Hopf bifurcation point and thus stabilises the focus and merges the stable node and the saddle node by a saddlenode bifurcation. Here, the Hop-bifurcation and the saddlenode bifurcation occur on different stationary state branches and at different noise levels. Moreover, the eigenfrequency ν of the oscillation on the lower branch decreases monotonically with increasing noise level. These properties explain the simulation results presented in Figure2(A): the γ-activity observed at low noise levels represents a nonlinear limit cycle above the Hopf-bifurcation point, while the low-amplitude γ-activity at larger noise levels represents a quasi-cycle about the stable focus below the Hopf-bifurcation point.

Figure 5

 5 Figure 5. Finite-size impact on oscillation frequency. The spectral power of the quasi-cycle on the lower branch shown in Fig.2(B)at D e /τ e = 0.4. The power distribution has been computed applying the Bartlett-Welch method for simulated signals V (t) = ∑ N i=1 V i (t)/N (black) as solution of Eq. (1) and x 1 (t) (red) as solution of Eqs.(9). Other parameters were identical to parameters used in Fig.2.

Figure 6 .

 6 Figure 6. Illustrative finite-size impact on transfer function and nonlinear gain. (A) Finite-size correction of the transfer function ∆G from Eq. (13) and its nonlinear gain ∆G from Eq.(14). For illustration reasons it is N = 1. (B) The transfer function and its nonlinear gain with fluctuations (N = 50, dotted line) and without fluctuations (∆G, ∆G = 0, solid line). Here, σ 2 = 0.1.

eV = x 1 φ e 1 = x 1 φ e 1 ) 2 .S 2 S 2

 11222 -(t-t )/τ e dη e k (t ) + ψ e † k Q k and Eq. (B8) reads t-t )/τ e ψ e, † k dY e (t ) + ψ e † k Q k t-t )/τ e dY e (t )w e (t) t-t )/τ e dR e (t ) ψ e † k Q k with dR e (t) = e dY e (t)/N. By the definition of Q k , we find in addition of the new matrices F , M have elements F i j = F 0 (X/cN -1/N), M i j = M 0 (X/cN -1/N) with the Bernoulli random number X. Consequently F , M have zero mean with variances identical to matrices F, M.In sum,V = (x 1 -Z(t))φ e 1 + w e (t) +F S 1 [ V] -M S 2 [ W]. (B9)Since we have assumed that V evolves slowly in time and does not fluctuate randomly, we may choose V = x 1 φ e 1 = x 1 e. The same line of argumentation holds for vector W yielding W = y 1 e with y 1 ∈ R. ThenS 1 [ V] = eS 1 [x 1 ], S 2 [ W] = eS 2 [y 1 ]. In addition (F e) i = U F , (M e) i = U M ∀i = 1, . . . , N with the normal distributed random variable U F ∼ N (0, σ 2 F /N), U M ∼ N (0, σ 2M /N) and we gain the final resultV = (x 1 -Z(t))e +w e (t) +U F eS 1 [x 1 ] -U M eS 2 [y 1 ] .(B10The terms Z(t), U F , U m emerge due to the network's finite size and vanish for N → ∞. In this case,V = u = xe + w and 1 + w n ) = S(x) + ξ (t)with the random number ξ (t) ∼ N (0, Var(S)/N). The last expression reflects the fundamental rule about the uncertainty of a sample mean. The theoretical mean value for an infinitely large sample is S(x) = lim + w)p(w)dw with the probability density function p(w) of the random numbers w n . Then Var(S(x)) = limN→∞ (x + w n ) -S2 (u) (x + w)p e (w)dw -S2 (x) . Since S(x) = S 0 Θ(x)is a Heaviside function, one can write S 2 (x) = S 2 0 S(x) and we finally gain Var(S(x)) = S 2 0 S(x) -S2 (x) and e N S[u] = S(x) + S 2 0 S(x) -S2 (x) N X with the random number X ∼ N (0, 1). Specifically S(x) = S 0 2 1 + erf x 2(D e /τ e )

ACKNOWLEDGMENTS

We thank the National Sciences and Engineering Research Council of Canada (NSERC Grants RGPIN-2017-06662 to JL), and INRIA (Action Exploratoire A/D Drugs to AH) for support of this research.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Eigenvalue spectrum of random Bernouilli matrices

A N × N random Bernoulli matrix K with entries K i j = X/cN, where the Bernoulli-distributed random numbers X have mean c and variance c(1c), has the mean value µ = 1/N and variance

Then the edge spectrum has the mean λ 1 = 1, which is also the maximum eigenvalue, and the variance

1 → 0 for infinite-size networks. Moreover, the bulk spectrum has zero mean and the bulk spectrum's maximum eigenvalue is λ k>1,max = 2 (1c)/cN. For infinite size networks, this maximum eigenvalue vanishes. Figure 7 shows the eigenvalues of such a matrix K pointing out both spectra. 1. For a random Bernoulli-matrix K ∈ R N×N as defined in section A with probability c and cN → ∞, the eigenvalue problem

In addition, |λ k>1 | → 0 for N → ∞.

For the bi-orthonormal basis {ψ

the corresponding eigenvalue problem reads

with ψ 1 = e/N. The symbol † indicates the transposition and complex-conjugate.

Proof.

1. According to [START_REF] Ding | Spectral distributions of adjacency and laplacian matrices of random graphs[END_REF][START_REF] Füredi | The eigenvalues of random symmetric matrices[END_REF][START_REF] Guilonnet | Bernoulli random matrices[END_REF][START_REF] May | Stability and complexity in model ecosystems[END_REF] , the spectrum of a random matrix includes a bulk spectrum and an edge spectrum. Typically the edge spectrum is localised far from the bulk spectrum, distributed around its mean value, and exists only if the expectation value of the random matrix does not vanish, e.g. see Fig. 2 in [START_REF] Hutt | Coherence resonance in random erdos-renyi neural networks : mean-field theory[END_REF] and Fig. 7. If the random matrix K has Bernoulli-distributed entries K i j = K 0 /cN • X with Bernoulli-distributed random numbers X with mean c and variance c(1c), then the mean value of K is µ = K 0 /N and its variance

We define that the first eigenvalues λ k of the edge spectrum with the first modes. Since there is a single eigenvalue in the edge spectrum, its mean is

and variance σ 2 1 . For symmetric Bernoulli matrices, σ 2 1 = 2σ 2 K . The edge spectrums eigenvectors have random entries with identity mean, i.e. φ 1 ≈ e. Since ∑ j=1 K i j (φ 1 ) j = λ 1 (φ 1 ) i ∀i = 1, . . . , N