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

Jérémie Lefebvre, Axel Hutt

► To cite this version:

Jérémie Lefebvre, Axel Hutt. Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2023, 10.1063/5.0167771. hal-04162581v2

HAL Id: hal-04162581

https://inria.hal.science/hal-04162581v2

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study

Jérémie Lefebvre, Axel Hutt

► To cite this version:

Jérémie Lefebvre, Axel Hutt. Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2023, 10.1063/5.0167771. hal-04162581v2

HAL Id: hal-04162581

https://inria.hal.science/hal-04162581v2

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study

J. Lefebvre<fnref n="1" />, <fnref n="2" />, <fnref n="3" /> and A. Hutt*<fnref n="4" />

Event-related synchronisation and desynchronisation (ERS/ERD) are well-known features found experimentally in brain signals during cognitive tasks. Their understanding promises to have much better insights into neural information processes in cognition. Under the hypothesis that neural information affects the endogenous neural noise level in populations, we propose to employ a stochastic mean-field model to explain ERS/ERD in the $\gamma$-frequency range. The work extends previous mean-field studies by deriving novel effects from finite network size. Moreover, numerical simulations of ERS/ERD and their analytical explanation by the mean-field model suggests several endogenous noise modulation schemes which may modulate the system's synchronisation.

Endogenous neuroelectric fluctuations are assumed to contribute or reflect information processing in the brain. We hypothesize that these fluctuations may induce event-related synchronisation and desynchronisation (ERS/ERD), which are spectral power enhancements and reductions, respectively, in certain frequency bands. A stochastic network model considers such endogenous noise and exhibits ERS/ERD in the $\gamma$-frequency range dependent on the endogenous noise level. Mathematical analysis and numerical simulations reveal the underlying mechanism called Additive-Noise Induced System Evolution (ANISE) in line with experimental results. The work proposes a novel perspective to explain ERS/ERD in various frequency bands to reveal underlying mechanisms of brain information processing.

## 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¹⁴), in time (e.g. coherent oscillations in given frequency bands¹⁶), or a combination of both (e.g. spatio-temporal dynamics such as travelling waves⁵⁰). 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 ¹⁹,²⁵,²⁹. This self-organisation is observed across a wide variety of complex systems²⁸,³⁸,³⁹,⁴⁶, 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²¹,³⁷, 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⁵², which determines their stability and resilience to perturbations³⁴.

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¹⁵. 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³⁰, in the corticothalamic feedback system³² or under various forms of neurostimulation⁴².

Additive input noise may further be used to account for the effect of stimuli and exogeneous inputs. For instance, high-frequency 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⁵⁸. 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⁴⁷. 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

<footnote n="1">¹) Krembil Brain Institute, University Health Network, Toronto, ON M5T 0S8, Canada</footnote>

<footnote n="2">²) Department of Biology, University of Ottawa, Ottawa, ON K1N 6N5, Canada</footnote>

<footnote n="3">³) Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada</footnote>

<footnote n="4">⁴) ICube, MLMS, University of Strasbourg; MIMESIS Team, Inria Nancy - Grand Est, 67000 Strasbourg, France</footnote>

<footnote>(*Electronic mail: axel.hutt@inria.fr)</footnote>

<footnote>(Dated: February 23, 2024)</footnote>
# Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study

J. Lefebvre<fnref n="1" />, <fnref n="2" />, <fnref n="3" /> and A. Hutt*<fnref n="4" />

Event-related synchronisation and desynchronisation (ERS/ERD) are well-known features found experimentally in brain signals during cognitive tasks. Their understanding promises to have much better insights into neural information processes in cognition. Under the hypothesis that neural information affects the endogenous neural noise level in populations, we propose to employ a stochastic mean-field model to explain ERS/ERD in the $\gamma$-frequency range. The work extends previous mean-field studies by deriving novel effects from finite network size. Moreover, numerical simulations of ERS/ERD and their analytical explanation by the mean-field model suggests several endogenous noise modulation schemes which may modulate the system's synchronisation.

Endogenous neuroelectric fluctuations are assumed to contribute or reflect information processing in the brain. We hypothesize that these fluctuations may induce event-related synchronisation and desynchronisation (ERS/ERD), which are spectral power enhancements and reductions, respectively, in certain frequency bands. A stochastic network model considers such endogenous noise and exhibits ERS/ERD in the $\gamma$-frequency range dependent on the endogenous noise level. Mathematical analysis and numerical simulations reveal the underlying mechanism called Additive-Noise Induced System Evolution (ANISE) in line with experimental results. The work proposes a novel perspective to explain ERS/ERD in various frequency bands to reveal underlying mechanisms of brain information processing.

## 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¹⁴), in time (e.g. coherent oscillations in given frequency bands¹⁶), or a combination of both (e.g. spatio-temporal dynamics such as travelling waves⁵⁰). 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 ¹⁹,²⁵,²⁹. This self-organisation is observed across a wide variety of complex systems²⁸,³⁸,³⁹,⁴⁶, 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²¹,³⁷, 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⁵², which determines their stability and resilience to perturbations³⁴.

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¹⁵. 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³⁰, in the corticothalamic feedback system³² or under various forms of neurostimulation⁴².

Additive input noise may further be used to account for the effect of stimuli and exogeneous inputs. For instance, high-frequency 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⁵⁸. 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⁴⁷. 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

<footnote n="1">¹) Krembil Brain Institute, University Health Network, Toronto, ON M5T 0S8, Canada</footnote>

<footnote n="2">²) Department of Biology, University of Ottawa, Ottawa, ON K1N 6N5, Canada</footnote>

<footnote n="3">³) Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada</footnote>

<footnote n="4">⁴) ICube, MLMS, University of Strasbourg; MIMESIS Team, Inria Nancy - Grand Est, 67000 Strasbourg, France</footnote>

<footnote>(*Electronic mail: axel.hutt@inria.fr)</footnote>

<footnote>(Dated: February 23, 2024)</footnote>
α-frequency band (8 – 12Hz), ERD reflects a state of cortical activation or enhanced neuronal processing associated with the cognitive or motor task⁴⁴. 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⁴⁰. 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³⁶ or during neurostimulation⁶³. 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³²,⁵⁸ and thus resembles coherence resonance³⁵,⁴⁹. 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²⁰,²²,⁵⁵,⁶⁰. Conversely, there is strong evidence that γ-rhythms originate from local circuits⁵¹,⁶¹. 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³⁰,⁵⁴. 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¹⁸,³⁶,⁴⁷,⁵⁷. To reduce the complexity of corresponding random network models, mean-field 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²⁶,³⁰,³¹,³⁵,⁵³ have considered primarily the case of infinite size, whereas finite-size networks may introduce additional random effects³¹,³⁵ 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 mean-field 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⁴² or coherence resonance³,³⁵. 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

<figure><img src="image_2.png" /><figcaption>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.⁵⁶ licensed under a Creative Commons Attribution (CC BY) license.</figcaption></figure>

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 \in \mathbb{R}^N$ and $W \in \mathbb{R}^N$, respectively, which represent mean dendritic currents of the corresponding micro-column (also known as minicolumn⁵,⁸,⁴⁵). The associated dynamics of these micro-columns obey

$$ 
\begin{aligned}
\tau_e dV &= (-V + FS_1[V] - MS_2[W] + eI_e)dt + dY_e(t) \\
\tau_i dW &= (-W + MS_1[V] - FS_2[W] + eI_i)dt + dY_i(t)
\end{aligned}
\quad (1) $$

with the N-dimensional identity vector $e = (1, \dots, 1)^T$, the intra-subnetwork random connectivity matrix $F$, the inter-subnetwork random connectivity matrix $M$ and the synaptic time scales $\tau_e$ and $\tau_i$. The connectivity matrices $F$ and $M$ are $N \times N$-dimensional random matrices. Both matrices share the structure of a random matrix $K \in \mathbb{R}^{N \times N}$, whose elements $\{(K)_{ij}\}$ are Bernoulli-distributed with probability $c$ and are equal to $(K)_{ij} = XK_o/cN$, $K_o > 0$, where $X$ is a Bernoulli random variable. Then the mean connectivities are $\mu_F = F_0/N$, $\mu_M = M_0/N$ with respective variance $\sigma_F^2 = F_o^2(1-c)/cN^2$, $\sigma_M^2 = M_o^2(1-c)/cN^2$.

The external stochastic input is a N-dimensional Wiener process $Y_{e,i} \in \mathbb{R}^N$, for which

$$ 
\begin{aligned}
\langle (dY_{e,i})_n(t) \rangle &= 0 \\
\langle (dY_{e,i})_n(t) (dY_{e,i})_m(t') \rangle &= 2D_{e,i} \delta_{nm} \delta(t-t') dt dt'
\end{aligned}
$$
α-frequency band (8 – 12Hz), ERD reflects a state of cortical activation or enhanced neuronal processing associated with the cognitive or motor task⁴⁴. 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⁴⁰. 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³⁶ or during neurostimulation⁶³. 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³²,⁵⁸ and thus resembles coherence resonance³⁵,⁴⁹. 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²⁰,²²,⁵⁵,⁶⁰. Conversely, there is strong evidence that γ-rhythms originate from local circuits⁵¹,⁶¹. 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³⁰,⁵⁴. 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¹⁸,³⁶,⁴⁷,⁵⁷. To reduce the complexity of corresponding random network models, mean-field 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²⁶,³⁰,³¹,³⁵,⁵³ have considered primarily the case of infinite size, whereas finite-size networks may introduce additional random effects³¹,³⁵ 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 mean-field 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⁴² or coherence resonance³,³⁵. 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

<figure><img src="image_2.png" /><figcaption>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.⁵⁶ licensed under a Creative Commons Attribution (CC BY) license.</figcaption></figure>

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 \in \mathbb{R}^N$ and $W \in \mathbb{R}^N$, respectively, which represent mean dendritic currents of the corresponding micro-column (also known as minicolumn⁵,⁸,⁴⁵). The associated dynamics of these micro-columns obey

$$ 
\begin{aligned}
\tau_e dV &= (-V + FS_1[V] - MS_2[W] + eI_e)dt + dY_e(t) \\
\tau_i dW &= (-W + MS_1[V] - FS_2[W] + eI_i)dt + dY_i(t)
\end{aligned}
\quad (1)
$$

with the N-dimensional identity vector $e = (1, \dots, 1)^T$, the intra-subnetwork random connectivity matrix $F$, the inter-subnetwork random connectivity matrix $M$ and the synaptic time scales $\tau_e$ and $\tau_i$. The connectivity matrices $F$ and $M$ are $N \times N$-dimensional random matrices. Both matrices share the structure of a random matrix $K \in \mathbb{R}^{N \times N}$, whose elements $\{(K)_{ij}\}$ are Bernoulli-distributed with probability $c$ and are equal to $(K)_{ij} = XK_o/cN$, $K_o > 0$, where $X$ is a Bernoulli random variable. Then the mean connectivities are $\mu_F = F_0/N$, $\mu_M = M_0/N$ with respective variance $\sigma_F^2 = F_o^2(1-c)/cN^2$, $\sigma_M^2 = M_o^2(1-c)/cN^2$.

The external stochastic input is a N-dimensional Wiener process $Y_{e,i} \in \mathbb{R}^N$, for which

$$ 
\begin{aligned}
\langle (dY_{e,i})_n(t) \rangle &= 0 \\
\langle (dY_{e,i})_n(t) (dY_{e,i})_m(t') \rangle &= 2D_{e,i} \delta_{nm} \delta(t-t') dt dt'
\end{aligned}
$$
for $n,m = 1, \dots, N$ and the noise variance $D_{e,i}$. The symbol $\langle \cdot \rangle$ denotes the ensemble average, $\delta_{nm}$ is the Kronecker symbol and $\delta(\cdot)$ denotes the $\delta$-distribution.

Nodes in the network exhibit a nonlinear transfer function $S_1(\cdot)$ and $S_2(\cdot)$ 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$^{1,27}$ with $S_1(x) = S_0\Theta(x)$, $S_0 > 1$ and $S_2(x) = \Theta(x)$ with the Heaviside function $\Theta(x) = 0 \forall x < 0$, $\Theta(x) = 1 \forall x \ge 0$. The function vector $S_{1,2}[x] \in \mathbb{R}^N$ is defined 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) = \sum_{n=1}^{N} x_n(t) \phi_n^e \quad (2)$$

$$W(t) = \sum_{n=1}^{N} y_n(t) \phi_n^i$$

with the vector basis $\{\phi_n^e\}$, $\{\phi_n^i\}$ and $\phi_n^e, \phi_n^i \in \mathbb{C}^N$ and a biorthonormal basis $\{\psi_k^e\}$, $\{\psi_k^i\}$ and $\psi_k^e, \psi_k^i \in \mathbb{C}^N$ with

$$\begin{aligned} \psi_k^{e, \dagger} \phi_l^e &= \delta_{kl}, & \psi_k^{i, \dagger} \phi_l^i &= \delta_{kl}, \\ \psi_k^{x, \dagger} \phi_l^y &= 0 \quad \forall x, y = e, i, \end{aligned}$$

Proposition 1 stipulates

$$\phi_1^e = \phi_1^i = \mathbf{e}, \quad \psi_1^e = \psi_1^i = \mathbf{e}/N.$$

Moreover, we define the eigenvalue problems

$$\psi_k^{e,i,\dagger} \mathbf{F} = \lambda_k^F \psi_k^{e,i,\dagger}, \quad \psi_k^{e,i,\dagger} \mathbf{M} = \lambda_k^M \psi_k^{e,i,\dagger} \quad (3)$$

for $k = 1, \dots, N$.

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

$$\begin{aligned} \tau_e dx_1 &= \left(-x_1 + \frac{\lambda_1^F}{N} e^T \mathbf{S}_1[\mathbf{V}] - \frac{\lambda_1^M}{N} e^T \mathbf{S}_2[\mathbf{W}] + I_e\right) dt + d\eta^e(t) \\ \tau_i dy_1 &= \left(-y_1 + \frac{\lambda_1^M}{N} e^T \mathbf{S}_1[\mathbf{V}] - \frac{\lambda_1^F}{N} e^T \mathbf{S}_2[\mathbf{W}] + I_i\right) dt + d\eta^i(t) \end{aligned} \quad (4)$$

with $d\eta^{e,i}(t) = e^T d\mathbf{Y}_{e,i}(t)/N$ and the remaining modes obey

$$\begin{aligned} \tau_e dx_k(t) &= \left(-x_k(t) + \lambda_k^F \psi_k^{e,\dagger} \mathbf{S}_1[\mathbf{V}] - \lambda_k^M \psi_k^{e,\dagger} \mathbf{S}_2[\mathbf{W}] + \psi_k^{e,\dagger} \mathbf{e} I_e\right) dt \\ &\quad + d\eta_k^e(t) \\ \tau_i dy_k(t) &= \left(-y_k(t) + \lambda_k^M \psi_k^{i,\dagger} \mathbf{S}_1[\mathbf{V}] - \lambda_k^F \psi_k^{i,\dagger} \mathbf{S}_2[\mathbf{W}] + \psi_k^{i,\dagger} \mathbf{e} I_i\right) dt \\ &\quad + d\eta_k^i(t) \end{aligned} \quad (5)$$

for $k > 1$ with $d\eta_k^{e,i}(t) = \psi_k^{e,i,\dagger} d\mathbf{Y}_{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

$$\begin{aligned} w_e(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{Y}_e(t') \\ Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{R}^e(t') \end{aligned}$$

and the finite size random corrections $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$. The finite size-fluctuations $R^e(t)$ are Wiener processes with $dR^e \sim \mathcal{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[\mathbf{u}] = (S_1(u_1), \dots, S_1(u_N))^\top \in \mathbb{R}^N$ like in Eq. (1) and $\mathbf{u} = x\mathbf{e} + \mathbf{w} \in \mathbb{R}^N$ with $x \in \mathbb{R}$ and random vector $\mathbf{w} \in \mathbb{R}^N$, $w_i \sim \mathcal{N}(0, D_e/\tau_e)$, Proposition 2 yields

$$\frac{\mathbf{e}^\top}{N} S[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} \chi(t)$$

with random variable $\chi(t) \sim \mathcal{N}(0, 1)$ and

$$\begin{aligned} \bar{S}(x) &= \int_{-\infty}^{\infty} S_1(x+w) p_e(w) dw \\ &= \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right) . \end{aligned}$$

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., $\phi_1^e \approx \mathbf{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 $\phi_1^e \neq \mathbf{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.
for $n,m = 1, \dots, N$ and the noise variance $D_{e,i}$. The symbol $\langle \cdot \rangle$ denotes the ensemble average, $\delta_{nm}$ is the Kronecker symbol and $\delta(\cdot)$ denotes the $\delta$-distribution.

Nodes in the network exhibit a nonlinear transfer function $S_1(\cdot)$ and $S_2(\cdot)$ 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$^{1,27}$ with $S_1(x) = S_0\Theta(x)$, $S_0 > 1$ and $S_2(x) = \Theta(x)$ with the Heaviside function $\Theta(x) = 0 \forall x < 0$, $\Theta(x) = 1 \forall x \ge 0$. The function vector $S_{1,2}[x] \in \mathbb{R}^N$ is defined 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 $\mathbf{F}$ and $\mathbf{M}$, will prove helpful. Assuming the basis composition

$$ \mathbf{V}(t) = \sum_{n=1}^{N} x_n(t) \phi_n^e \quad (2) $$

$$ \mathbf{W}(t) = \sum_{n=1}^{N} y_n(t) \phi_n^i $$

with the vector basis $\{\phi_n^e\}$, $\{\phi_n^i\}$ and $\phi_n^e, \phi_n^i \in \mathbb{C}^N$ and a biorthonormal basis $\{\psi_k^e\}$, $\{\psi_k^i\}$ and $\psi_k^e, \psi_k^i \in \mathbb{C}^N$ with

$$ \begin{aligned} \psi_k^{e, \dagger} \phi_l^e &= \delta_{kl}, & \psi_k^{i, \dagger} \phi_l^i &= \delta_{kl}, \\ \psi_k^{x, \dagger} \phi_l^y &= 0 \quad \forall x, y = e, i, \end{aligned} $$

Proposition 1 stipulates

$$ \phi_1^e = \phi_1^i = \mathbf{e}, \quad \psi_1^e = \psi_1^i = \mathbf{e}/N. $$

Moreover, we define the eigenvalue problems

$$ \psi_k^{e,i,\dagger} \mathbf{F} = \lambda_k^F \psi_k^{e,i,\dagger}, \quad \psi_k^{e,i,\dagger} \mathbf{M} = \lambda_k^M \psi_k^{e,i,\dagger} \quad (3) $$

for $k = 1, \dots, N$.

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

$$ \begin{aligned} \tau_e dx_1 &= \left(-x_1 + \frac{\lambda_1^F}{N} e^\mathrm{T} \mathbf{S}_1[\mathbf{V}] - \frac{\lambda_1^M}{N} e^\mathrm{T} \mathbf{S}_2[\mathbf{W}] + I_e\right) dt + d\eta^e(t) \\ \tau_i dy_1 &= \left(-y_1 + \frac{\lambda_1^M}{N} e^\mathrm{T} \mathbf{S}_1[\mathbf{V}] - \frac{\lambda_1^F}{N} e^\mathrm{T} \mathbf{S}_2[\mathbf{W}] + I_i\right) dt + d\eta^i(t) \end{aligned} \quad (4) $$

with $d\eta^{e,i}(t) = e^\mathrm{T} d\mathbf{Y}_{e,i}(t)/N$ and the remaining modes obey

$$ \begin{aligned} \tau_e dx_k(t) &= \left(-x_k(t) + \lambda_k^F \psi_k^{e,\dagger} \mathbf{S}_1[\mathbf{V}] - \lambda_k^M \psi_k^{e,\dagger} \mathbf{S}_2[\mathbf{W}] + \psi_k^{e,\dagger} e I_e\right) dt \\ &\quad + d\eta_k^e(t) \\ \tau_i dy_k(t) &= \left(-y_k(t) + \lambda_k^M \psi_k^{i,\dagger} \mathbf{S}_1[\mathbf{V}] - \lambda_k^F \psi_k^{i,\dagger} \mathbf{S}_2[\mathbf{W}] + \psi_k^{i,\dagger} e I_i\right) dt \\ &\quad + d\eta_k^i(t) \end{aligned} \quad (5) $$

for $k > 1$ with $d\eta_k^{e,i}(t) = \psi_k^{e,i,\dagger} d\mathbf{Y}_{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 $\mathbf{V}(t)$ in Eq. (2) as

$$ \mathbf{V} = (x_1 - Z(t) + U_F S_1[x_1] - U_M S_2[y_1])e + \mathbf{w}_e(t) $$

with

$$ \begin{aligned} \mathbf{w}_e(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{Y}_e(t') \\ Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{R}^e(t') \end{aligned} $$

and the finite size random corrections $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$. The finite size-fluctuations $R^e(t)$ are Wiener processes with $dR^e \sim \mathcal{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 $\mathbf{S}[\mathbf{u}] = (S_1(u_1), \dots, S_1(u_N))^\mathrm{T} \in \mathbb{R}^N$ like in Eq. (1) and $\mathbf{u} = x\mathbf{e} + \mathbf{w} \in \mathbb{R}^N$ with $x \in \mathbb{R}$ and random vector $\mathbf{w} \in \mathbb{R}^N$, $w_i \sim \mathcal{N}(0, D_e/\tau_e)$, Proposition 2 yields

$$ \frac{\mathbf{e}^\mathrm{T}}{N} \mathbf{S}[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} \chi(t) $$

with random variable $\chi(t) \sim \mathcal{N}(0, 1)$ and

$$ \begin{aligned} \bar{S}(x) &= \int_{-\infty}^{\infty} S_1(x+w) p_e(w) dw \\ &= \frac{S_0}{2} \left( 1 + \erf \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right). \end{aligned} $$

The dynamics of $\mathbf{W}(t)$ can also be written similarly and the derivation resembles closely to the one of $\mathbf{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., $\phi_1^e \approx \mathbf{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 $\phi_1^e \neq \mathbf{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.
<figure><img src="image_3.png" /><figcaption>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\tau_e$ and changed to $D_e = 0.8\tau_e$ at $t = 5$s and changed back to $D_e = 0.25\tau_e$ at $t = 15$s. Bottom panel: Time-frequency plot of instantaneous spectral power (windowed Fourier Transform convolved with a Hamming window with 2s window width and overlap of 1.8s). (B) Top panel: the stationary state $\bar{x}_1$ dependent on $D_e/\tau_e$ revealing s stable node (dotted line), saddle node (dashed line) and a stable focus (solid line). Bottom panel: the maximum eigenvalue real part $\lambda_m$ (left hand side) and the eigenfrequency of the focus (right hand side). Other parameters were $D_i/\tau_i = 0.2$, $N = 100$, $F_0 = 2.17$, $M_0 = 3.87$, $\tau_e = 0.005$, $\tau_i = 0.02$, $S_0 = 1.7$. and the integration time step was 0.5ms utilising the Euler-Maruyama method.</figcaption></figure>

The results of Proposition 2 allow us to recast Eqs. (4) as

$$ 
\begin{aligned}
\tau_e dx_1 &= (-x_1 + \lambda_1^F G_e (x_1 - Z(t) + T_1(x_1, y_1), t) \\
&\quad - \lambda_1^M G_i (y_1 - W(t) + T_2(x_1, y_1), t) + I_e) dt + dR^e(t) \\
\tau_i dy_1 &= (-y_1 + \lambda_1^M G_e (x_1 - Z(t) + T_1(x_1, y_1), t) \\
&\quad - \lambda_1^F G_i (y_1 - W(t) + T_2(x_1, y_1), t) + I_i) dt + dR^i(t)
\end{aligned}
\qquad (6) $$

with

$$ 
\begin{aligned}
G_e(x_1, t) &= \bar{G}_e(x_1) + \Delta G_e(x_1) \cdot X_e(t) \\
G_i(y_1, t) &= \bar{G}_i(y_1) + \Delta G_i(y_1) \cdot X_i(t),
\end{aligned}
\qquad (7) $$

for temporal fluctuations $X_{e,i}(t) \sim \mathcal{N}(0, 1)$ and where

$$ 
\begin{aligned}
\bar{G}_e(x_1) &= \int_{-\infty}^{\infty} S_1(x_1 + w) p_e(w) dw \\
\bar{G}_i(y_1) &= \int_{-\infty}^{\infty} S_2(y_1 + w) p_i(w) dw \\
\Delta G_e(x_1) &= \sqrt{\frac{S_0^2 \bar{G}_e(x_1) - \bar{G}_e^2(x_1)}{N}} \\
\Delta G_i(y_1) &= \sqrt{\frac{\bar{G}_i(y_1) - \bar{G}_i^2(y_1)}{N}}.
\end{aligned}
$$

Combining previous results, the network exhibits finite-size

fluctuations

$$ 
\begin{aligned}
dR^e &= \xi_e dt, & \xi_e &\sim \mathcal{N}(0, D_e/N) \\
dR^i &= \xi_i dt, & \xi_i &\sim \mathcal{N}(0, D_i/N) \\
Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} dR^e(t') \\
W(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} dR^i(t'),
\end{aligned}
$$

where random state-dependent corrections of the nonlinear interactions occur through the terms $\Delta G_{e,i}$, which now possess random state-dependent thresholds

$$ 
\begin{aligned}
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),
\end{aligned}
$$

with $U_F \sim \mathcal{N}(0, F_0^2(1-c)/cN^3)$, $U_M \sim \mathcal{N}(0, M_0^2(1-c)/cN^3)$. For large networks, $F_0^2(1-c)/cN^3 \ll D_{e,i}/N$, $M_0^2(1-c)/cN^3 \ll D_{e,i}/N$ and thus the state-dependent thresholds are negligible compared to the finite-size fluctuations. If $N \to \infty$, $T_{1,2} \to 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) = \frac{1}{\sqrt{2\pi D_{e,i} / \tau_{e,i}}} e^{-w^2 / 2(D_{e,i} / \tau_{e,i})} $$
<figure><img src="image_3.png" /><figcaption>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\tau_e$ and changed to $D_e = 0.8\tau_e$ at $t = 5s$ and changed back to $D_e = 0.25\tau_e$ at $t = 15s$. Bottom panel: Time-frequency plot of instantaneous spectral power (windowed Fourier Transform convolved with a Hamming window with 2s window width and overlap of 1.8s). (B) Top panel: the stationary state $\bar{x}_1$ dependent on $D_e/\tau_e$ revealing s stable node (dotted line), saddle node (dashed line) and a stable focus (solid line). Bottom panel: the maximum eigenvalue real part $\lambda_m$ (left hand side) and the eigenfrequency of the focus (right hand side). Other parameters were $D_i/\tau_i = 0.2$, $N = 100$, $F_0 = 2.17$, $M_0 = 3.87$, $\tau_e = 0.005$, $\tau_i = 0.02$, $S_0 = 1.7$. and the integration time step was 0.5ms utilising the Euler-Maruyama method.</figcaption></figure>

The results of Proposition 2 allow us to recast Eqs. (4) as

$$ 
\begin{aligned}
\tau_e dx_1 &= (-x_1 + \lambda_1^F G_e (x_1 - Z(t) + T_1(x_1, y_1), t) \\
&\quad - \lambda_1^M G_i (y_1 - W(t) + T_2(x_1, y_1), t) + I_e) dt + dR^e(t) \\
\tau_i dy_1 &= (-y_1 + \lambda_1^M G_e (x_1 - Z(t) + T_1(x_1, y_1), t) \\
&\quad - \lambda_1^F G_i (y_1 - W(t) + T_2(x_1, y_1), t) + I_i) dt + dR^i(t)
\end{aligned}
\quad (6) $$

with

$$ 
\begin{aligned}
G_e(x_1, t) &= \bar{G}_e(x_1) + \Delta G_e(x_1) \cdot X_e(t) \\
G_i(y_1, t) &= \bar{G}_i(y_1) + \Delta G_i(y_1) \cdot X_i(t),
\end{aligned}
\quad (7) $$

for temporal fluctuations $X_{e,i}(t) \sim \mathcal{N}(0, 1)$ and where

$$ 
\begin{aligned}
\bar{G}_e(x_1) &= \int_{-\infty}^{\infty} S_1(x_1 + w) p_e(w) dw \\
\bar{G}_i(y_1) &= \int_{-\infty}^{\infty} S_2(y_1 + w) p_i(w) dw \\
\Delta G_e(x_1) &= \sqrt{\frac{S_0^2 \bar{G}_e(x_1) - \bar{G}_e^2(x_1)}{N}} \\
\Delta G_i(y_1) &= \sqrt{\frac{\bar{G}_i(y_1) - \bar{G}_i^2(y_1)}{N}}.
\end{aligned}
$$

Combining previous results, the network exhibits finite-size

fluctuations

$$ 
\begin{aligned}
dR^e &= \xi_e dt, \quad \xi_e \sim \mathcal{N}(0, D_e/N) \\
dR^i &= \xi_i dt, \quad \xi_i \sim \mathcal{N}(0, D_i/N) \\
Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} dR^e(t') \\
W(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} dR^i(t'),
\end{aligned}
$$

where random state-dependent corrections of the nonlinear interactions occur through the terms $\Delta G_{e,i}$, which now possess random state-dependent thresholds

$$ 
\begin{aligned}
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),
\end{aligned}
$$

with $U_F \sim \mathcal{N}(0, F_0^2(1-c)/cN^3)$, $U_M \sim \mathcal{N}(0, M_0^2(1-c)/cN^3)$. For large networks, $F_0^2(1-c)/cN^3 \ll D_{e,i}/N$, $M_0^2(1-c)/cN^3 \ll D_{e,i}/N$ and thus the state-dependent thresholds are negligible compared to the finite-size fluctuations. If $N \to \infty$, $T_{1,2} \to 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) = \frac{1}{\sqrt{2\pi D_{e,i} / \tau_{e,i}}} e^{-w^2 / 2(D_{e,i} / \tau_{e,i})} $$
and consequently, according to Proposition 2,

$$ 
\begin{aligned}
\bar{G}_e(x) &= \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right) \\
\bar{G}_i(x) &= \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_i/\tau_i)}} \right) \right).
\end{aligned}
\quad (8) $$

### III. SIMULATIONS

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

$$ 
\begin{aligned}
dx_1 &= \left( -\frac{x_1}{\tau_e} + \frac{\lambda_1^F}{\tau_e} G_e(x_1 - Z(t), t) - \frac{\lambda_1^M}{\tau_e} G_i(y_1 - W(t), t) + \frac{J_1^e}{\tau_e} \right) dt + \frac{1}{\tau_e} dR^e(t) \\
dy_1 &= \left( -\frac{y_1}{\tau_i} + \frac{\lambda_1^M}{\tau_i} G_e(x_1 - Z(t), t) - \frac{\lambda_1^F}{\tau_i} G_i(y_1 - W(t), t) + \frac{J_1^i}{\tau_i} \right) dt + \frac{1}{\tau_i} dR^i(t)
\end{aligned}
\quad (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 $\sqrt{D_e}$ and $\sqrt{D_i}$, respectively, while the nonlinear gain of the transfer functions $G_e$ and $G_i$ (their slope) decrease nonlinearly with $\sqrt{D_e}$ and $\sqrt{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$^{29,33}$.

#### A. Synchronisation modulated by noise

Experimental studies$^{18,36,57}$ 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.

<figure><img src="image_4.png" /><figcaption>Figure 3. γ-desynchronisation by decreasing excitatory input noise (A) Top panel: the time-dependent change of the noise variance $D_e/\tau_e$ in the simulation. Bottom panel: $x_1(t)$ simulated from Eq. (9) dependent on $D_e/\tau_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/\tau_e$ moves the system's stationary state $\bar{x}_1$ 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/\tau_e$ reflecting a hysteresis. The period on the top branch reflects a non-oscillatory, 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.</figcaption></figure>

#### 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 time-frequency plot of the numerical mean $\bar{V}(t) = \sum_{i=1}^N 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 $\bar{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 $\bar{x}_1$ resulting from the fixed point condition $dx_1/dt = dy_1/dt = 0$ in the absence of finite-size noise, i.e. $N \to \infty$. Linearizing Eqs. (9) about the stationary states $\bar{x}_1, \bar{y}_1$ with $x_1(t) = \bar{x}_1 + u(t)$, $y_1 = \bar{y}_1 + v(t)$ we obtain the following dynamical system

$$ \tau dz/dt = Lz + s(t) \quad (10) $$

with $z = (u, v)^T$ and

$$ \tau = \begin{pmatrix} \tau_e & 0 \\ 0 & \tau_i \end{pmatrix}, \quad L = \begin{pmatrix} \lambda_1^F G'_e - 1 & -\lambda_1^M G'_i \\ \lambda_1^M G'_e & -\lambda_1^F G'_i - 1 \end{pmatrix} \quad (11) $$

$$ s = \begin{pmatrix} -\lambda_1^F G'_e Z(t) + \lambda_1^M G'_i W(t) \\ -\lambda_1^M G'_e Z(t) + \lambda_1^F G'_i W(t) \end{pmatrix}. \quad (12) $$

Here, $G'_{e,i} = dG_{e,i}(x)/dx$ is evaluated at $x = \bar{x}_1, \bar{y}_1$. The stationary states' stability is determined by the maximum eigenvalue real part $\lambda_m$ of $L$. If the eigenvalue $\lambda_m$ is complex-valued with $\mathcal{Im}(\lambda_m) = 2\pi\nu$, then $\nu$ is called eigenfrequency.
and consequently, according to Proposition 2,

$$ \bar{G}_e(x) = \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right) \quad (8) $$

$$ \bar{G}_i(x) = \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_i/\tau_i)}} \right) \right). $$

### III. SIMULATIONS

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

$$ \begin{aligned} dx_1 &= \left( -\frac{x_1}{\tau_e} + \frac{\lambda_1^F}{\tau_e} G_e(x_1 - Z(t), t) - \frac{\lambda_1^M}{\tau_e} G_i(y_1 - W(t), t) \right. \\ & \qquad \left. + \frac{J_1^e}{\tau_e} \right) dt + \frac{1}{\tau_e} dR^e(t) \\ dy_1 &= \left( -\frac{y_1}{\tau_i} + \frac{\lambda_1^M}{\tau_i} G_e(x_1 - Z(t), t) - \frac{\lambda_1^F}{\tau_i} G_i(y_1 - W(t), t) \right. \\ & \qquad \left. + \frac{J_1^i}{\tau_i} \right) dt + \frac{1}{\tau_i} dR^i(t) \end{aligned} \quad (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 $\sqrt{D_e}$ and $\sqrt{D_i}$, respectively, while the nonlinear gain of the transfer functions $G_e$ and $G_i$ (their slope) decrease nonlinearly with $\sqrt{D_e}$ and $\sqrt{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$^{29,33}$.

#### A. Synchronisation modulated by noise

Experimental studies$^{18,36,57}$ 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.

<figure><img src="image_4.png" /><figcaption>Figure 3. γ-desynchronisation by decreasing excitatory input noise (A) Top panel: the time-dependent change of the noise variance $D_e/\tau_e$ in the simulation. Bottom panel: $x_1(t)$ simulated from Eq. (9) dependent on $D_e/\tau_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/\tau_e$ moves the system's stationary state $\bar{x}_1$ 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/\tau_e$ reflecting a hysteresis. The period on the top branch reflects a non-oscillatory, 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.</figcaption></figure>

#### 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 time-frequency plot of the numerical mean $\bar{V}(t) = \sum_{i=1}^N 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 $\bar{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 $\bar{x}_1$ resulting from the fixed point condition $dx_1/dt = dy_1/dt = 0$ in the absence of finite-size noise, i.e. $N \to \infty$. Linearizing Eqs. (9) about the stationary states $\bar{x}_1, \bar{y}_1$ with $x_1(t) = \bar{x}_1 + u(t)$, $y_1 = \bar{y}_1 + v(t)$ we obtain the following dynamical system

$$ \tau dz/dt = Lz + s(t) \quad (10) $$

with $z = (u, v)^T$ and

$$ \tau = \begin{pmatrix} \tau_e & 0 \\ 0 & \tau_i \end{pmatrix}, \quad L = \begin{pmatrix} \lambda_1^F G'_e - 1 & -\lambda_1^M G'_i \\ \lambda_1^M G'_e & -\lambda_1^F G'_i - 1 \end{pmatrix} \quad (11) $$

$$ s = \begin{pmatrix} -\lambda_1^F G'_e Z(t) + \lambda_1^M G'_i W(t) \\ -\lambda_1^M G'_e Z(t) + \lambda_1^F G'_i W(t) \end{pmatrix}. \quad (12) $$

Here, $G'_{e,i} = dG_{e,i}(x)/dx$ is evaluated at $x = \bar{x}_1, \bar{y}_1$. The stationary states' stability is determined by the maximum eigenvalue real part $\lambda_m$ of $L$. If the eigenvalue $\lambda_m$ is complex-valued with $\mathcal{Im}(\lambda_m) = 2\pi\nu$, then $\nu$ is called eigenfrequency.
<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

Figure 2(B,top panel) shows the stationary state $\bar{x}_1$ 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 saddle-node bifurcation. Here, the Hop-furcation and the saddle-node bifurcation occur on different stationary state branches and at different noise levels. Moreover, the eigenfrequency $v$ 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

<figure><img src="image_6.png" /><figcaption>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 $\bar{V}(t) = \sum_{i=1}^{N} 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.</figcaption></figure>

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.

### ERD 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⁴⁷. 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-
<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

Figure 2(B,top panel) shows the stationary state $\bar{x}_1$ 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 $\lambda_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 saddle-node bifurcation. Here, the Hop-furcation and the saddle-node bifurcation occur on different stationary state branches and at different noise levels. Moreover, the eigenfrequency $v$ 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

<figure><img src="image_6.png" /><figcaption>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/\tau_e = 0.4$. The power distribution has been computed applying the Bartlett-Welch method for simulated signals $\bar{V}(t) = \sum_{i=1}^{N} 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.</figcaption></figure>

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.

### ERD 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⁴⁷. 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-
<figure><img src="image_7.png" /><figcaption>Figure 6. Illustrative finite-size impact on transfer function and nonlinear gain. (A) Finite-size correction of the transfer function $\Delta G$ from Eq. (13) and its nonlinear gain $\Delta 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 ($\Delta G, \Delta G' = 0$, solid line). Here, $\sigma^2 = 0.1$.</figcaption></figure>

quotencies compared to large networks. Hence the network size determines the system's principal frequencies and their distribution.

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

$$ 
\begin{aligned}
G(V) &= \bar{G}(V) + \Delta G(V) \cdot X \\
\bar{G}(V) &= \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{V}{\sqrt{2}\sigma} \right) \right) \\
\Delta G(V) &= \sqrt{\frac{\bar{G}(V) - \bar{G}^2(V)}{N}}
\end{aligned} 
\quad (13) $$

with the random variable $X \sim \mathcal{N}(0, 1)$ and real parameter $\sigma > 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) = \bar{G}'(V) \left( 1 + \underbrace{\frac{1 - 2\bar{G}(V)}{2\sqrt{\bar{G}(V) - \bar{G}^2(V)N}}}_{=\Delta G'(V)} \cdot X \right) . \quad (14) $$

by taking the derivative of $G(V)$ in Eq. (13). Figure 6(A) presents the finite-size correction to the transfer function $\Delta G$ from Eq. (13) and its derivative $\Delta G'$ from Eq. (14). The absolute value of the nonlinear gain finite-size correction is maximum about $V = \pm 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 $\gamma$-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 = \pm 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 mean-field 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 non-coherent behavior subjected to the endogenous noise level and thus exhibits synchronisation and desynchronisation, respectively (Figs. 2-4). Conversely, smaller networks induce finite-size 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 $\gamma$-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$^{2,7}$ 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$^{6,12}$. Moreover, previous other work has provided evidence that ERD may result from reduced network inhibition$^{59}$. 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.
<figure><img src="image_7.png" /><figcaption>Figure 6. Illustrative finite-size impact on transfer function and nonlinear gain. (A) Finite-size correction of the transfer function $\Delta G$ from Eq. (13) and its nonlinear gain $\Delta 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 ($\Delta G, \Delta G' = 0$, solid line). Here, $\sigma^2 = 0.1$.</figcaption></figure>

quotencies compared to large networks. Hence the network size determines the system's principal frequencies and their distribution.

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

$$ 
\begin{aligned}
G(V) &= \bar{G}(V) + \Delta G(V) \cdot X \\
\bar{G}(V) &= \frac{1}{2} \left( 1 + \operatorname{erf} \left( \frac{V}{\sqrt{2}\sigma} \right) \right) \\
\Delta G(V) &= \sqrt{\frac{\bar{G}(V) - \bar{G}^2(V)}{N}}
\end{aligned}
\quad (13) $$

with the random variable $X \sim \mathcal{N}(0, 1)$ and real parameter $\sigma > 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) = \bar{G}'(V) \left( 1 + \underbrace{\frac{1 - 2\bar{G}(V)}{2\sqrt{\bar{G}(V) - \bar{G}^2(V)N}}}_{=\Delta G'(V)} \cdot X \right) . \quad (14) $$

by taking the derivative of $G(V)$ in Eq. (13). Figure 6(A) presents the finite-size correction to the transfer function $\Delta G$ from Eq. (13) and its derivative $\Delta G'$ from Eq. (14). The absolute value of the nonlinear gain finite-size correction is maximum about $V = \pm 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 $\gamma$-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 = \pm 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 mean-field 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 non-coherent behavior subjected to the endogenous noise level and thus exhibits synchronisation and desynchronisation, respectively (Figs. 2-4). Conversely, smaller networks induce finite-size 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 $\gamma$-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$^{2,7}$ 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$^{6,12}$. Moreover, previous other work has provided evidence that ERD may result from reduced network inhibition$^{59}$. 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<sup>23,41</sup> or small-world networks<sup>4</sup> 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<sup>9,24</sup>.

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<sup>10,64</sup>. 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<sup>48</sup>.

## DATA AVAILABILITY STATEMENT

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

## 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.

## Appendix A: Eigenvalue spectrum of random Bernoulli matrices

A $N \times N$ random Bernoulli matrix $\mathbf{K}$ with entries $K_{ij} = X/cN$, where the Bernoulli-distributed random numbers $X$ have mean $c$ and variance $c(1-c)$, has the mean value $\mu = 1/N$ and variance $\sigma_K^2 = (1-c)/cN^2$. Then the edge spectrum has the mean $\lambda_1 = 1$, which is also the maximum eigenvalue, and the variance $\sigma_1^2 = 2(1-c)/cN^2$. Clearly, $\lambda_1 = 1$, $\sigma_1^2 \to 0$ for infinite-size networks. Moreover, the bulk spectrum has zero mean and the bulk spectrum's maximum eigenvalue is $\lambda_{k>1,max} = 2\sqrt{(1-c)/cN}$. For infinite size networks, this maximum eigenvalue vanishes. Figure 7 shows the eigenvalues of such a matrix $\mathbf{K}$ pointing out both spectra.

## Appendix B: Propositions and their proofs

**Proposition 1.** The spectral properties of random matrices are considered in two parts.

<figure><img src="image_8.png" /><figcaption>Figure 7. Edge and bulk spectrum of a random Bernoulli matrix. $N = 100, c = 0.9$</figcaption></figure>

1. For a random Bernoulli-matrix $\mathbf{K} \in \mathbb{R}^{N \times N}$ as defined in section A with probability $c$ and $cN \to \infty$, the eigenvalue problem

$$ \mathbf{K}\phi_k = \lambda_k\phi_k, \quad k = 1, \dots, N $$

with complex eigenvalues $\lambda_k \in \mathbb{C}$ and right hand-side eigenvectors $\phi_k \in \mathbb{C}^N$ is solved by $\lambda_1 = K_0$ and $\phi_1 = \mathbf{e}_1$. In addition, $|\lambda_{k>1}| \to 0$ for $N \to \infty$.

2. For the bi-orthonormal basis $\{\psi_k\}$, $\psi_k \in \mathbb{C}^N$ with

$$ \psi_k^\dagger \phi_l = \delta_{kl}, \quad (B1) $$

the corresponding eigenvalue problem reads

$$ \psi_k^\dagger \mathbf{K} = \lambda_k \psi_k^\dagger, \quad k = 1, \dots, N. $$

with $\psi_1 = \mathbf{e}/N$. The symbol $\dagger$ indicates the transposition and complex-conjugate.

**Proof.** 1. According to<sup>11,13,17,43</sup>, 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<sup>35</sup> and Fig. 7. If the random matrix $\mathbf{K}$ has Bernoulli-distributed entries $K_{ij} = K_0/cN \cdot X$ with Bernoulli-distributed random numbers $X$ with mean $c$ and variance $c(1-c)$, then the mean value of $\mathbf{K}$ is $\mu = K_0/N$ and its variance $\sigma_K^2 = K_0^2 c(1-c)/c^2 N^2 = K_0^2 (1-c)/cN^2$. We define that the first eigenvalues $\lambda_k$ of the edge spectrum with the first modes. Since there is a single eigenvalue in the edge spectrum, its mean is

$$ \lambda_1 = N\mu = K_0 \quad (B2) $$

and variance $\sigma_1^2$. For symmetric Bernoulli matrices, $\sigma_1^2 = 2\sigma_K^2$. The edge spectrums eigenvectors have random entries with identity mean, i.e. $\phi_1 \approx \mathbf{e}$. Since

$$ \sum_{j=1} K_{ij}(\phi_1)_j = \lambda_1(\phi_1)_i \quad \forall i = 1, \dots, N $$
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<sup>23,41</sup> or small-world networks<sup>4</sup> 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<sup>9,24</sup>.

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<sup>10,64</sup>. 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<sup>48</sup>.

## DATA AVAILABILITY STATEMENT

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

## 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.

## Appendix A: Eigenvalue spectrum of random Bernoulli matrices

A $N \times N$ random Bernoulli matrix $\mathbf{K}$ with entries $K_{ij} = X/cN$, where the Bernoulli-distributed random numbers $X$ have mean $c$ and variance $c(1-c)$, has the mean value $\mu = 1/N$ and variance $\sigma_K^2 = (1-c)/cN^2$. Then the edge spectrum has the mean $\lambda_1 = 1$, which is also the maximum eigenvalue, and the variance $\sigma_1^2 = 2(1-c)/cN^2$. Clearly, $\lambda_1 = 1$, $\sigma_1^2 \to 0$ for infinite-size networks. Moreover, the bulk spectrum has zero mean and the bulk spectrum's maximum eigenvalue is $\lambda_{k>1,max} = 2\sqrt{(1-c)/cN}$. For infinite size networks, this maximum eigenvalue vanishes. Figure 7 shows the eigenvalues of such a matrix $\mathbf{K}$ pointing out both spectra.

## Appendix B: Propositions and their proofs

**Proposition 1.** The spectral properties of random matrices are considered in two parts.

<figure><img src="image_8.png" /><figcaption>Figure 7. Edge and bulk spectrum of a random Bernoulli matrix. $N = 100, c = 0.9$</figcaption></figure>

1. For a random Bernoulli-matrix $\mathbf{K} \in \mathbb{R}^{N \times N}$ as defined in section A with probability $c$ and $cN \to \infty$, the eigenvalue problem

$$ \mathbf{K}\phi_k = \lambda_k\phi_k, \quad k = 1, \dots, N $$

with complex eigenvalues $\lambda_k \in \mathbb{C}$ and right hand-side eigenvectors $\phi_k \in \mathbb{C}^N$ is solved by $\lambda_1 = K_0$ and $\phi_1 = \mathbf{e}_1$. In addition, $|\lambda_{k>1}| \to 0$ for $N \to \infty$.

2. For the bi-orthonormal basis $\{\psi_k\}$, $\psi_k \in \mathbb{C}^N$ with

$$ \psi_k^\dagger \phi_l = \delta_{kl}, \quad (B1) $$

the corresponding eigenvalue problem reads

$$ \psi_k^\dagger \mathbf{K} = \lambda_k \psi_k^\dagger, \quad k = 1, \dots, N. $$

with $\psi_1 = \mathbf{e}/N$. The symbol $\dagger$ indicates the transposition and complex-conjugate.

**Proof.** 1. According to<sup>11,13,17,43</sup>, 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<sup>35</sup> and Fig. 7. If the random matrix $\mathbf{K}$ has Bernoulli-distributed entries $K_{ij} = K_0/cN \cdot X$ with Bernoulli-distributed random numbers $X$ with mean $c$ and variance $c(1-c)$, then the mean value of $\mathbf{K}$ is $\mu = K_0/N$ and its variance $\sigma_K^2 = K_0^2 c(1-c)/c^2 N^2 = K_0^2 (1-c)/cN^2$. We define that the first eigenvalues $\lambda_k$ of the edge spectrum with the first modes. Since there is a single eigenvalue in the edge spectrum, its mean is

$$ \lambda_1 = N\mu = K_0 \quad (B2) $$

and variance $\sigma_1^2$. For symmetric Bernoulli matrices, $\sigma_1^2 = 2\sigma_K^2$. The edge spectrums eigenvectors have random entries with identity mean, i.e. $\phi_1 \approx \mathbf{e}$. Since

$$ \sum_{j=1} K_{ij}(\phi_1)_j = \lambda_1(\phi_1)_i \quad \forall i = 1, \dots, N $$
and $\sigma_K^2, \sigma_I^2 \to 0$ for $N \to \infty$, it follows that $(\phi_1)_i \to 1$. This last statement is based on the corresponding results for symmetric Bernoulli matrices for very large networks with $cN \to \infty$¹³ only, but is confirmed by numerical verification.

The mean eigenvalue of the bulk spectrum is $E[\lambda_k] = 0$ with maximum value$^{43,62}$

$$ \lambda_{k>1,\max} = 2\sigma_K \sqrt{N} \quad (B3) $$

$$ = 2K_0 \sqrt{\frac{1-c}{cN}} \quad (B4) $$

and the entries of their eigenvectors are random$^{11,17}$. We observe that $|\lambda_{k>1,\max}| \to 0$ for $N \to \infty$.

2. For the orthonormal basis {$\psi$}, we assume

$$ \psi_k^\dagger \mathbf{K} = \lambda_k' \psi_k^\dagger $$

and we obtain

$$ \begin{aligned} \mathbf{K}\phi_k &= \lambda_k\phi_k \\ \psi_k^\dagger \mathbf{K}\phi_k &= \lambda_k \psi_k^\dagger \phi_k \\ \psi_k^\dagger \mathbf{K}\phi_k &= \lambda_k \\ \lambda_k' \psi_k^\dagger \phi_k &= \lambda_k \\ \lambda_k' &= \lambda_k. \end{aligned} $$

Since $\phi_1 = \mathbf{e}$, we obtain $\psi_1 = \mathbf{e}/N$ by relation (B1). □

<figure><img src="image_9.png" /></figure>

**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

$$ \mathbf{V} = (x_1 - Z(t) + U_F S_1[x_1] - U_M S_2[y_1])e + \mathbf{w}_e(t) \quad (B5) $$

with

$$ \begin{aligned} \mathbf{w}_e(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{Y}_e(t') \\ Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{R}^e(t'). \end{aligned} $$

and the finite size random corrections $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$. The finite size-fluctuations $\mathbf{R}^e(t)$ are Wiener processes with $d\mathbf{R}^e \sim \mathcal{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 $\mathbf{S}[\mathbf{u}] = (S_1(u_1), \dots, S_1(u_N))^\top \in \mathbb{R}^N$ like in Eq. (1) and $\mathbf{u} = \mathbf{x}\mathbf{e} + \mathbf{w} \in \mathbb{R}^N$ with $\mathbf{x} \in \mathcal{R}$ and random vector $\mathbf{w} \in \mathcal{R}^N$, $w_i \sim \mathcal{N}(0, D_e/\tau_e)$. Then

$$ \frac{\mathbf{e}^\top}{N} \mathbf{S}[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} \chi(t) \quad (B6) $$

with random variable $\chi(t) \sim \mathcal{N}(0, 1)$ and

$$ \begin{aligned} \bar{S}(x) &= \int_{-\infty}^{\infty} S_1(x+w) p_e(w) dw \\ &= \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right). \end{aligned} $$

Proof. 1. We consider

$$ \tau_e dx_k(t) = (-x_k(t) + \psi_k^{e\dagger} \mathbf{Q}_k) dt + d\eta_k^e(t) \quad (B7) $$

for $k = 2, \dots, N$ motivated by Eq. (5) with

$$ \mathbf{Q}_k(t) = \lambda_k^F \mathbf{S}_1[\mathbf{V}(t)] - \lambda_k^M \mathbf{S}_2[\mathbf{W}(t)] + \mathbf{e}I_e. $$

The solution of (B7) is

$$ x_k(t) = \int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\eta_k^e(t') + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} \mathbf{Q}_k(t') dt' $$

and Eq. (2) reads then

$$ \mathbf{V} = x_1 \phi_1^e + \sum_{k=2}^{N} \phi_k^e \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} d\mathbf{Y}_e(t') \\ + \sum_{k=2}^{N} \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} \mathbf{Q}_k(t') dt'. \quad (B8) $$

This equation is implicit in $\mathbf{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 $\mathbf{V}(t)$ evolves very slowly over time, i.e. slower than the time scale $\tau_e$ and thus $\mathbf{Q}(t) \approx \text{const}$. This implies that $\mathbf{Q}(t)$ does not fluctuate randomly. Consequently,

$$ \mathbf{Q}_k = \lambda_k^F \mathbf{S}_1[\bar{\mathbf{V}}] - \lambda_k^M \mathbf{S}_2[\bar{\mathbf{W}}] + \mathbf{e}I_e $$

and $\mathbf{V} \approx \bar{\mathbf{V}}$, $\mathbf{W} \approx \bar{\mathbf{W}}$ is constant in time. Then the mode projection amplitudes $x_k(t)$ obey Ornstein-Uhlenbeck processes

$$ x_k(t) = \int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\eta_k^e(t') + \psi_k^{e\dagger} \mathbf{Q}_k $$
and $\sigma_K^2, \sigma_I^2 \to 0$ for $N \to \infty$, it follows that $(\phi_1)_i \to 1$. This last statement is based on the corresponding results for symmetric Bernoulli matrices for very large networks with $cN \to \infty$¹³ only, but is confirmed by numerical verification.

The mean eigenvalue of the bulk spectrum is $E[\lambda_k] = 0$ with maximum value$^{43,62}$

$$ \lambda_{k>1,\max} = 2\sigma_K \sqrt{N} \quad (B3) $$

$$ = 2K_0 \sqrt{\frac{1-c}{cN}} \quad (B4) $$

and the entries of their eigenvectors are random$^{11,17}$. We observe that $|\lambda_{k>1,\max}| \to 0$ for $N \to \infty$.

2. For the orthonormal basis {$\psi$}, we assume

$$ \psi_k^\dagger \mathbf{K} = \lambda_k' \psi_k^\dagger $$

and we obtain

$$ \begin{aligned} \mathbf{K}\phi_k &= \lambda_k\phi_k \\ \psi_k^\dagger \mathbf{K}\phi_k &= \lambda_k \psi_k^\dagger \phi_k \\ \psi_k^\dagger \mathbf{K}\phi_k &= \lambda_k \\ \lambda_k' \psi_k^\dagger \phi_k &= \lambda_k \\ \lambda_k' &= \lambda_k. \end{aligned} $$

Since $\phi_1 = \mathbf{e}$, we obtain $\psi_1 = \mathbf{e}/N$ by relation (B1). $\square$

<figure><img src="image_9.png" /></figure>

**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

$$ \mathbf{V} = (x_1 - Z(t) + U_F S_1[x_1] - U_M S_2[y_1])e + \mathbf{w}_e(t) \quad (B5) $$

with

$$ \begin{aligned} \mathbf{w}_e(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{Y}_e(t') \\ Z(t) &= \frac{1}{\tau} \int_{-\infty}^{t} e^{-(t-t')/\tau} d\mathbf{R}^e(t'). \end{aligned} $$

and the finite size random corrections $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$. The finite size-fluctuations $\mathbf{R}^e(t)$ are Wiener processes with $d\mathbf{R}^e \sim \mathcal{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 $\mathbf{S}[\mathbf{u}] = (S_1(u_1), \dots, S_1(u_N))^\top \in \mathbb{R}^N$ like in Eq. (1) and $\mathbf{u} = \mathbf{x}\mathbf{e} + \mathbf{w} \in \mathbb{R}^N$ with $\mathbf{x} \in \mathbb{R}$ and random vector $\mathbf{w} \in \mathbb{R}^N$, $w_i \sim \mathcal{N}(0, D_e/\tau_e)$. Then

$$ \frac{\mathbf{e}^\top}{N} \mathbf{S}[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} \chi(t) \quad (B6) $$

with random variable $\chi(t) \sim \mathcal{N}(0, 1)$ and

$$ \begin{aligned} \bar{S}(x) &= \int_{-\infty}^{\infty} S_1(x+w) p_e(w) dw \\ &= \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right). \end{aligned} $$

Proof. 1. We consider

$$ \tau_e dx_k(t) = (-x_k(t) + \psi_k^{e\dagger} \mathbf{Q}_k) dt + d\eta_k^e(t) \quad (B7) $$

for $k = 2, \dots, N$ motivated by Eq. (5) with

$$ \mathbf{Q}_k(t) = \lambda_k^F \mathbf{S}_1[\mathbf{V}(t)] - \lambda_k^M \mathbf{S}_2[\mathbf{W}(t)] + \mathbf{e}I_e. $$

The solution of (B7) is

$$ x_k(t) = \int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\eta_k^e(t') + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} \mathbf{Q}_k(t') dt' $$

and Eq. (2) reads then

$$ \mathbf{V} = x_1 \phi_1^e + \sum_{k=2}^{N} \phi_k^e \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} d\mathbf{Y}_e(t') + \sum_{k=2}^{N} \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e\dagger} \mathbf{Q}_k(t') dt'. \quad (B8) $$

This equation is implicit in $\mathbf{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 $\mathbf{V}(t)$ evolves very slowly over time, i.e. slower than the time scale $\tau_e$ and thus $\mathbf{Q}(t) \approx \text{const}$. This implies that $\mathbf{Q}(t)$ does not fluctuate randomly. Consequently,

$$ \mathbf{Q}_k = \lambda_k^F \mathbf{S}_1[\bar{\mathbf{V}}] - \lambda_k^M \mathbf{S}_2[\bar{\mathbf{W}}] + \mathbf{e}I_e $$

and $\mathbf{V} \approx \bar{\mathbf{V}}$, $\mathbf{W} \approx \bar{\mathbf{W}}$ is constant in time. Then the mode projection amplitudes $x_k(t)$ obey Ornstein-Uhlenbeck processes

$$ x_k(t) = \int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\eta_k^e(t') + \psi_k^{e\dagger} \mathbf{Q}_k $$
and Eq. (B8) reads

$$
\begin{align*}
\mathbf{V} &= x_1 \phi_1^e + \sum_{k=2}^{N} \phi_k^e \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e,\dagger} d\mathbf{Y}_e(t') + \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \right) d\mathbf{Y}_e(t') + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \left( \mathbf{I} - \phi_1^e \psi_1^{e,\dagger} \right) d\mathbf{Y}_e(t') + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \underbrace{\int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\mathbf{Y}_e(t')}_{\mathbf{w}_e(t)} \\
&\quad - \phi_1^e \underbrace{\int_{-\infty}^{t} e^{-(t-t')/\tau_e} dR^e(t')}_{Z(t)} \\
&\quad + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k
\end{align*}
$$

with $dR^e(t) = \mathbf{e}^\mathrm{T} d\mathbf{Y}_e(t)/N$. By the definition of $\mathbf{Q}_k$, we
find in addition

$$
\begin{align*}
\sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \mathbf{Q}_k &= \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \lambda_k^F \right) \mathbf{S}_1[\bar{\mathbf{V}}] \\
&\quad - \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \lambda_k^M \right) \mathbf{S}_2[\bar{\mathbf{W}}] \\
&\quad + \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \right) \mathbf{e}I_e \\
&= \underbrace{\left( \mathbf{F} - \lambda_1^F \phi_1^e \psi_1^{e\dagger} \right)}_{\mathbf{F}'} \mathbf{S}_1[\bar{\mathbf{V}}] \\
&\quad - \underbrace{\left( \mathbf{M} - \lambda_1^M \phi_1^e \psi_1^{e\dagger} \right)}_{\mathbf{M}'} \mathbf{S}_2[\bar{\mathbf{W}}] \\
&\quad + \underbrace{\left( \mathbf{I} - \phi_1^e \psi_1^{e\dagger} \right)}_{=0} \mathbf{e}I_e.
\end{align*}
$$

Here, we have used the identities

$$
\mathbf{F} = \sum_{k=1}^{N} \lambda_k^F \phi_k^e \psi_k^{e\dagger} , \quad \mathbf{M} = \sum_{k=1}^{N} \lambda_k^M \phi_k^e \psi_k^{e\dagger}
$$

and the elements of the new matrices **F', M'** have ele-
ments $F'_{ij} = F_0(X/cN - 1/N)$, $M'_{ij} = M_0(X/cN - 1/N)$
with the Bernoulli random number $X$. Consequently
**F', M'** have zero mean with variances identical to ma-
trices **F, M**.

In sum,

$$
\mathbf{V} = (x_1 - Z(t))\phi_1^e + \mathbf{w}_e(t) \\
+ \mathbf{F}'\mathbf{S}_1[\bar{\mathbf{V}}] - \mathbf{M}'\mathbf{S}_2[\bar{\mathbf{W}}]. \quad (B9)
$$

Since we have assumed that $\bar{\mathbf{V}}$ evolves slowly in time and does not fluctuate randomly, we may choose $\bar{\mathbf{V}} = x_1 \phi_1^e = x_1 \mathbf{e}$. The same line of argumentation holds for vector $\mathbf{W}$ yielding $\bar{\mathbf{W}} = y_1 \mathbf{e}$ with $y_1 \in \mathcal{R}$. Then $\mathbf{S}_1[\bar{\mathbf{V}}] = \mathbf{e} \mathbf{S}_1[y_1]$, $\mathbf{S}_2[\bar{\mathbf{W}}] = \mathbf{e} \mathbf{S}_2[y_1]$. In addition $(\mathbf{F}'\mathbf{e})_i = U_F$, $(\mathbf{M}'\mathbf{e})_i = U_M \quad \forall i = 1, \dots, N$ with the normal distributed random variable $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$ and we gain the final result

$$
\mathbf{V} = (x_1 - Z(t))e \\
+ \mathbf{w}_e(t) + U_F \mathbf{e} \mathbf{S}_1[x_1] - U_M \mathbf{e} \mathbf{S}_2[y_1]. \quad \text{(B10)}
$$

2. The terms $Z(t)$, $U_F$, $U_M$ emerge due to the network's finite size and vanish for $N \to \infty$. In this case, $\mathbf{V} = \mathbf{u} = \mathbf{x}\mathbf{e} + \mathbf{w}$ and

$$
\begin{align*}
\frac{\mathbf{e}^\mathrm{T}}{N} \mathbf{S}[\mathbf{u}] &= \frac{1}{N} \sum_{n=1}^{N} S(u_n) \\
&= \frac{1}{N} \sum_{n=1}^{N} S(x_1 + w_n) \\
&= \bar{S}(x) + \xi(t)
\end{align*}
$$

with the random number $\xi(t) \sim \mathcal{N}(0, \text{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

$$
\begin{align*}
\bar{S}(x) &= \lim_{N \to \infty} \frac{1}{N} \sum_{n=1}^{N} S(x + w_n) \\
&= \int_{-\infty}^{\infty} S(x + w) p(w) dw
\end{align*}
$$

with the probability density function $p(w)$ of the ran-
dom numbers $w_n$. Then

$$
\begin{align*}
\mathrm{Var}(S(x)) &= \lim_{N \to \infty} \frac{1}{N} \sum_{n=1}^{N} (S(x+w_n) - \bar{S}(u))^2 \\
&= \lim_{N \to \infty} \left( \frac{1}{N} \sum_{n=1}^{N} S^2(x+w_n) - \bar{S}^2(u) \right) \\
&= \int_{-\infty}^{\infty} S^2(x+w)p_e(w)dw - \bar{S}^2(x).
\end{align*}
$$

Since $S(x) = S_0\Theta(x)$ is a Heaviside function, one can
write $S^2(x) = S_0^2 S(x)$ and we finally gain

$$
\mathrm{Var}(S(x)) = S_0^2 \bar{S}(x) - \bar{S}^2(x)
$$

and

$$
\frac{\mathbf{e}^\mathrm{T}}{N} \mathbf{S}[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} X
$$

with the random number $X \sim \mathcal{N}(0, 1)$. Specifically

$$
\bar{S}(x) = \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right)
$$
and Eq. (B8) reads

$$
\begin{align*}
\mathbf{V} &= x_1 \phi_1^e + \sum_{k=2}^{N} \phi_k^e \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \psi_k^{e,\dagger} d\mathbf{Y}_e(t') + \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \right) d\mathbf{Y}_e(t') + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \int_{-\infty}^{t} e^{-(t-t')/\tau_e} \left( \mathbf{I} - \phi_1^e \psi_1^{e,\dagger} \right) d\mathbf{Y}_e(t') + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k \\
&= x_1 \phi_1^e + \underbrace{\int_{-\infty}^{t} e^{-(t-t')/\tau_e} d\mathbf{Y}_e(t')}_{\mathbf{w}_e(t)} \\
&\quad - \phi_1^e \underbrace{\int_{-\infty}^{t} e^{-(t-t')/\tau_e} dR^e(t')}_{Z(t)} \\
&\quad + \sum_{k=2}^{N} \phi_k^e \psi_k^{e,\dagger} \mathbf{Q}_k
\end{align*}
$$

with $dR^e(t) = \mathbf{e}^\mathrm{T} d\mathbf{Y}_e(t)/N$. By the definition of $\mathbf{Q}_k$, we find in addition

$$
\begin{align*}
\sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \mathbf{Q}_k &= \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \lambda_k^F \right) \mathbf{S}_1[\bar{\mathbf{V}}] \\
&\quad - \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \lambda_k^M \right) \mathbf{S}_2[\bar{\mathbf{W}}] \\
&\quad + \left( \sum_{k=2}^{N} \phi_k^e \psi_k^{e\dagger} \right) \mathbf{e}I_e \\
&= \underbrace{\left( \mathbf{F} - \lambda_1^F \phi_1^e \psi_1^{e\dagger} \right)}_{\mathbf{F}'} \mathbf{S}_1[\bar{\mathbf{V}}] \\
&\quad - \underbrace{\left( \mathbf{M} - \lambda_1^M \phi_1^e \psi_1^{e\dagger} \right)}_{\mathbf{M}'} \mathbf{S}_2[\bar{\mathbf{W}}] \\
&\quad + \underbrace{\left( \mathbf{I} - \phi_1^e \psi_1^{e\dagger} \right)}_{=0} \mathbf{e}I_e.
\end{align*}
$$

Here, we have used the identities

$$
\mathbf{F} = \sum_{k=1}^{N} \lambda_k^F \phi_k^e \psi_k^{e\dagger} , \quad \mathbf{M} = \sum_{k=1}^{N} \lambda_k^M \phi_k^e \psi_k^{e\dagger}
$$

and the elements of the new matrices **F', M'** have ele-
ments $F'_{ij} = F_0(X/cN - 1/N)$, $M'_{ij} = M_0(X/cN - 1/N)$
with the Bernoulli random number $X$. Consequently
**F', M'** have zero mean with variances identical to ma-
trices **F, M**.

In sum,

$$
\mathbf{V} = (x_1 - Z(t))\phi_1^e + \mathbf{w}_e(t) \\
+ \mathbf{F}'\mathbf{S}_1[\bar{\mathbf{V}}] - \mathbf{M}'\mathbf{S}_2[\bar{\mathbf{W}}]. \quad (\text{B9})
$$

Since we have assumed that $\bar{\mathbf{V}}$ evolves slowly in time and does not fluctuate randomly, we may choose $\bar{\mathbf{V}} = x_1 \phi_1^e = x_1 \mathbf{e}$. The same line of argumentation holds for vector $\mathbf{W}$ yielding $\bar{\mathbf{W}} = y_1 \mathbf{e}$ with $y_1 \in \mathcal{R}$. Then $\mathbf{S}_1[\bar{\mathbf{V}}] = \mathbf{e} \mathbf{S}_1[y_1]$, $\mathbf{S}_2[\bar{\mathbf{W}}] = \mathbf{e} \mathbf{S}_2[y_1]$. In addition $(\mathbf{F}'\mathbf{e})_i = U_F$, $(\mathbf{M}'\mathbf{e})_i = U_M \quad \forall i = 1, \dots, N$ with the normal distributed random variable $U_F \sim \mathcal{N}(0, \sigma_F^2/N)$, $U_M \sim \mathcal{N}(0, \sigma_M^2/N)$ and we gain the final result

$$
\mathbf{V} = (x_1 - Z(t))e \\
+ \mathbf{w}_e(t) + U_F \mathbf{e} \mathbf{S}_1[x_1] - U_M \mathbf{e} \mathbf{S}_2[y_1]. \quad \text{(B10)}
$$

2. The terms $Z(t)$, $U_F$, $U_M$ emerge due to the network's finite size and vanish for $N \to \infty$. In this case, $\mathbf{V} = \mathbf{u} = \mathbf{x}\mathbf{e} + \mathbf{w}$ and

$$
\begin{align*}
\frac{\mathbf{e}^\mathrm{T}}{N} \mathbf{S}[\mathbf{u}] &= \frac{1}{N} \sum_{n=1}^{N} S(u_n) \\
&= \frac{1}{N} \sum_{n=1}^{N} S(x_1 + w_n) \\
&= \bar{S}(x) + \xi(t)
\end{align*}
$$

with the random number $\xi(t) \sim \mathcal{N}(0, \text{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

$$
\begin{align*}
\bar{S}(x) &= \lim_{N \to \infty} \frac{1}{N} \sum_{n=1}^{N} S(x + w_n) \\
&= \int_{-\infty}^{\infty} S(x + w) p(w) dw
\end{align*}
$$

with the probability density function $p(w)$ of the ran-
dom numbers $w_n$. Then

$$
\begin{align*}
\mathrm{Var}(S(x)) &= \lim_{N \to \infty} \frac{1}{N} \sum_{n=1}^{N} (S(x+w_n) - \bar{S}(u))^2 \\
&= \lim_{N \to \infty} \left( \frac{1}{N} \sum_{n=1}^{N} S^2(x+w_n) - \bar{S}^2(u) \right) \\
&= \int_{-\infty}^{\infty} S^2(x+w)p_e(w)dw - \bar{S}^2(x).
\end{align*}
$$

Since $S(x) = S_0\Theta(x)$ is a Heaviside function, one can
write $S^2(x) = S_0^2 S(x)$ and we finally gain

$$
\mathrm{Var}(S(x)) = S_0^2 \bar{S}(x) - \bar{S}^2(x)
$$

and

$$
\frac{\mathbf{e}^{\mathrm{T}}}{N} \mathbf{S}[\mathbf{u}] = \bar{S}(x) + \sqrt{\frac{S_0^2 \bar{S}(x) - \bar{S}^2(x)}{N}} X
$$

with the random number $X \sim \mathcal{N}(0, 1)$. Specifically

$$
\bar{S}(x) = \frac{S_0}{2} \left( 1 + \operatorname{erf} \left( \frac{x}{\sqrt{2(D_e/\tau_e)}} \right) \right)
$$
with $w_i \sim \mathcal{N}(0, D_e/\tau_e)$.

<figure><img src="image_10.png" /></figure>
with $w_i \sim \mathcal{N}(0, D_e/\tau_e)$.

□
REFERENCES

¹M. Abeles. *Corticonics*. Cambridge University Press, 1991.

²C. Andrew and G. Pfurtscheller. Event-related coherence as a tool for studying dynamic interaction of brain regions. *Neurophysiol.*, 98:144–148, 1996.

³E. Baspinar, L. Schüler, S. Olmi, and A. Zakharova. Coherence resonance in neuronal populations: mean-field versus network model. *Phys. Rev. E*, 103(3-1):032308, 2021.

⁴D.S. Bassett and E.T. Bullmore. Small-world brain networks revisited. *Neuroscientist*, 23(5):499–516, 2017.

⁵M. Bennett. An attempt at a unified theory of the neocortical microcircuit in sensory cortex. *Front. Neural Circuits*, 14:40, 2020.

⁶M. Boly, E. Balteau, C. Schnakers, C. Degueldre G. Moonen, A. Luxen, C. Phillips P. Peigneux, P. Maquet, and S. Laureys. Baseline brain activity fluctuations predict somatosensory perception in humans. *Proc Natl Acad Sci USA*, 17:12187–12192, 2007.

⁷S.L. Bressler, R. Coppola, and R. Nakamura. Episodic multi-regional cortical coherence at multiple frequencies during visual task performance. *Nature*, 366:153–156, 1993.

⁸D.P. Buxhoeveden and M.F. Casanova. The minicolumn hypothesis in neuroscience. *Brain*, 125(5):935–951, 2002.

⁹A. Castelnovo, B. Graziano, F. Ferrarelli, and A. D'Agostino. Sleep spindles and slow waves in schizophrenia and related disorders: main findings, challenges and future perspectives. *Eur. J. Neurosci.*, 48:2738–2758, 2018.

¹⁰M.M. Churchland, B.M. Yu, J.P. Cunningham, L.P. Sugrue, M.R. Cohen, G.S. Corrado, W.T. Newsome, A.M. Clark, P. Hosseini, B.B. Scott, D.C. Bradley, M.A. Smith, A. Kohn, J.A. Movshon, K.M. Armstrong, T. Moore, S.W. Chang, L.H. Snyder, S.G. Lisberger, N.J. Priebe, I.M. Finn, D. Fester, S.I. Ryu, G. Santhanam, M. Sahani, and K.V. Shenoy. Stimulus onset quenches neural variability: a widespread cortical phenomenon. *Nat Neurosci.*, 13(3):369–378, 2010.

¹¹X. Ding and T. Jiang. Spectral distributions of adjacency and laplacian matrices of random graphs. *Ann. Appl. Prob.*, 20(6):2086–2117, 2010.

¹²M.D. Fox, A.Z. Snyder, J.L. Vincent, and M.E. Raichle. Intrinsic fluctuations within cortical systems account for intertrial variability in human behavior. *Neuron*, 56(1):171–184, 2007.

¹³Z. Füredi and J. Komlos. The eigenvalues of random symmetric matrices. *Combinatorica*, 1(3):233–241, 1981.

¹⁴R.J. Gardner, E. Hermansen, and M. Pachitariu. Toroidal topology of population activity in grid cells. *Nature*, 602:123–128, 2022.

¹⁵W. Gerstner and W. Kistler. *Spiking Neuron Models*. Cambridge University Press, Cambridge, 2002.

¹⁶J. González, M. Cavelli, A. Mondino, N. Rubido, A. Bl Tort, and P. Torterolo. Communication through coherence by means of cross-frequency coupling. *Neuroscience*, 449:157–164, 2020.

¹⁷A. Guilonnet. Bernoulli random matrices. submitted, 2021.

¹⁸A.R. Haig, E. Gordon, J.J. Wright, R.A. Meares, and H. Bahramali. Synchro nous cortical gamma-band activity in task-relevant cognition. *Neuroreport*, 11:669–675, 2000.

¹⁹H. Haken. *Brain Dynamics*. Springer, Berlin, 2002.

²⁰M. Halgren, I. Ulbert, H. Bastuji, D. Fabó, L. Erőss, M. Rey, O. Devinsky, W.K. Doyle, R. Mak-McCully, E. Halgren, L. Wittner, P. Chauvel, G. Heit, E. Eskandar, A. Mandell, and S.S. Cash. The generation and propagation of the human alpha rhythm. *Proc. Natl. Acad. Sci. USA*, 116(47):23772–23782, 2019.

²¹D. Hansel and H. Solominsky. Modeling feature selectivity in local cortical circuits. In C.Koch and I.Segev, editors, *Methods in Neuronal Modeling: From Ions to Networks*, chapter 13. MIT Press, Cambridge, 2 edition, 1998.

²²M. Hashemi, A. Hutt, and J. Sleigh. How the cortico-thalamic feedback affects the EEG power spectrum over frontal and occipital regions during propofol-induced anaesthetic sedation. *J. Comput. Neurosci.*, 39(1):155–159, 2015.

²³P. Holme. Rare and everywhere: Perspectives on scale-free networks. *Nat. Commun.*, 10:1016, 2019.

²⁴X. Huang, W.C. Troy, S.J. Schiff, Q. Yang, H. Ma, C.R. Laing, and J.Y. Wu. Spiral waves in disinhibited mammalian neocortex. *J.Neurosc.*, 24(44):9897–9902, 2004.

²⁵A. Hutt. Brain connectivity reduction reflects disturbed self-organisation of the brain: Neural disorders and general anaesthesia. In V. Cutsuridis, editor, *Multiscale models of Brain Disorders*, pages 207–218. Springer-Verlag, Berlin, 2019.

²⁶A. Hutt. Additive noise-induced system evolution (ANISE). *Front. Appl. Math. Stat.*, 8:879866, 2022.

²⁷A. Hutt and L. Buhry. Study of GABAergic extra-synaptic tonic inhibition in single neurons and neural populations by traversing neural scales: application to propofol-induced anaesthesia. *J. Comput. Neurosci.*, 37(3):417–437, 2014.

²⁸A. Hutt and H. Haken, editors. *Synergetics*. Springer-Verlag, New York, 2020.

²⁹A. Hutt and J. Lefebvre. Additive noise tunes the self-organization in complex systems. In A. Hutt and H. Haken, editors, *Synergetics*, Encyclopedia of Complexity and Systems Science Series, pages 183–196. Springer, New York, 2020.

³⁰A. Hutt and J. Lefebvre. Arousal fluctuations govern oscillatory transitions between dominant $\gamma$ and $\alpha$ occipital activity during eyes open/closed conditions. *Brain Topography*, 35:108–120, 2021.

³¹A. Hutt, J. Lefebvre, D. Hight, and H. Kaiser. Phase coherence induced by additive gaussian and non-gaussian noise in excitable networks with application to burst suppression-like brain signals. *Front. Appl. Math. Stat.*, 5:69, 2020.

³²A. Hutt, J. Lefebvre, D. Hight, and J. Sleigh. Suppression of underlying neuronal fluctuations mediates EEG slowing during general anaesthesia. *Neuroimage*, 179:414–428, 2018.

³³A. Hutt, A. Mierau, and J. Lefebvre. Dynamic control of synchronous activity in networks of spiking neurons. *PLoS One*, 11(9):e0161488, 2016.

³⁴A. Hutt, S. Rich, T.A. Valiante, and J. Lefebvre. Intrinsic neural diversity quenches the dynamic volatility of neural networks. *Proc. Nat. Acad. Sci. USA*, 120(28):e2218841120, 2023.

³⁵A. Hutt, T. Wahl, N. Voges, J. Hausmann, and J. Lefebvre. Coherence resonance in random erdós-renyi neural networks : mean-field theory. *Front. Appl. Math. Stat.*, 7:697904, 2021.

³⁶M. Iijima, R. Mase, M. Osawa, S. Shimizu, and S. Uchiyama. Event-related synchronization and desynchronization of high-frequency electroencephalographic activity during a visual Go/No-Go paradigm. *Neurophysiobiology*, 71(1):17–24, 2015.

³⁷M. Kaiser, C. C. Hilgetag, and R. Kötter. Hierarchy and dynamics of neural networks. *Front. Neuroinf.*, 4:112, 2010.

³⁸M. Kaschube, M. Schnabel, and F. Wolf. Self-organization and the selection of pinwheel density in visual cortical development. *New J. Phys.*, 10:015009, 2008.

³⁹J.A.S. Kelso. *Dynamic Patterns: The Self-Organization of Brain and Behavior*. MIT Press, Cambridge, 1995.

⁴⁰W. Klimesch, P. Sauseng, and S. Hanslmayr. EEG alpha oscillations: the inhibition-timing hypothesis. *Brain Res. Rev.*, 53(1):63–88, 2007.

⁴¹U. Lee, G. Oh, S. Kim, G. Noh, B. Choi, and G.A. Mashour. Brain networks maintain a scale-free organization across consciousness, anesthesia, and recovery: evidence for adaptive reconfiguration. *Anesthesiology*, 113(5):1081–1091, 2010.

⁴²J. Lefebvre, A. Hutt, and F. Frohlich. Stochastic resonance mediates the state-dependent effect of periodic stimulation on cortical alpha oscillations. *eLife*, 6:e32054, 2017.

⁴³R. M. May. *Stability and complexity in model ecosystems*. Princeton University Press, 2001.

⁴⁴A. Mierau, W. Klimesch, and J. Lefebvre. State-dependent alpha peak frequency shifts: Experimental evidence, potential mechanisms and functional implications. *Neuroscience*, 360:146–154, 2017.

⁴⁵V.B. Mountcastle. Modality and topographic properties of single neurons of cat's somatic sensory cortex. *Neurophysiol.*, 20:408

⁴⁶G. Nicolis and I. Prigogine. *Self-Organization in Non-Equilibrium Systems: From Dissipative Structures to Order Through Fluctuations*. J. Wiley and Sons, New York, 1977.

⁴⁷G. Pfurtscheller and F.H. Lopes da Silva. Event-related EEG/MEG synchronization and desynchronization: basic principles. *Clin. Neurophysiol.*, 110(11):1842–1857, 1999.

⁴⁸B. Pietras, V. Schmutz, and T. Schwalger. Mesoscopic description of hippocampal replay and metastability in spiking neural networks with short-term plasticity. *PLoS Comput Biol*, 18(12):e1010809, 2022.

⁴⁹A.S. Pikovsky and J. Kurths. Coherence resonance in a noise-driven excitable system. *Phys. Rev. Lett.*, 78:775–778, 1997.
REFERENCES

¹M. Abeles. *Corticonics*. Cambridge University Press, 1991.

²C. Andrew and G. Pfurtscheller. Event-related coherence as a tool for studying dynamic interaction of brain regions. *Neurophysiol.*, 98:144–148, 1996.

³E. Baspinar, L. Schüler, S. Olmi, and A. Zakharova. Coherence resonance in neuronal populations: mean-field versus network model. *Phys. Rev. E*, 103(3-1):032308, 2021.

⁴D.S. Bassett and E.T. Bullmore. Small-world brain networks revisited. *Neuroscientist*, 23(5):499–516, 2017.

⁵M. Bennett. An attempt at a unified theory of the neocortical microcircuit in sensory cortex. *Front. Neural Circuits*, 14:40, 2020.

⁶M. Boly, E. Balteau, C. Schnakers, C. Degueldre G. Moonen, A. Luxen, C. Phillips P. Peigneux, P. Maquet, and S. Laureys. Baseline brain activity fluctuations predict somatosensory perception in humans. *Proc Natl Acad Sci USA*, 17:12187–12192, 2007.

⁷S.L. Bressler, R. Coppola, and R. Nakamura. Episodic multi-regional cortical coherence at multiple frequencies during visual task performance. *Nature*, 366:153–156, 1993.

⁸D.P. Buxhoeveden and M.F. Casanova. The minicolumn hypothesis in neuroscience. *Brain*, 125(5):935–951, 2002.

⁹A. Castelnovo, B. Graziano, F. Ferrarelli, and A. D'Agostino. Sleep spindles and slow waves in schizophrenia and related disorders: main findings, challenges and future perspectives. *Eur. J. Neurosci.*, 48:2738–2758, 2018.

¹⁰M.M. Churchland, B.M. Yu, J.P. Cunningham, L.P. Sugrue, M.R. Cohen, G.S. Corrado, W.T. Newsome, A.M. Clark, P. Hosseini, B.B. Scott, D.C. Bradley, M.A. Smith, A. Kohn, J.A. Movshon, K.M. Armstrong, T. Moore, S.W. Chang, L.H. Snyder, S.G. Lisberger, N.J. Priebe, I.M. Finn, D. Fester, S.I. Ryu, G. Santhanam, M. Sahani, and K.V. Shenoy. Stimulus onset quenches neural variability: a widespread cortical phenomenon. *Nat Neurosci.*, 13(3):369–378, 2010.

¹¹X. Ding and T. Jiang. Spectral distributions of adjacency and laplacian matrices of random graphs. *Ann. Appl. Prob.*, 20(6):2086–2117, 2010.

¹²M.D. Fox, A.Z. Snyder, J.L. Vincent, and M.E. Raichle. Intrinsic fluctuations within cortical systems account for intertrial variability in human behavior. *Neuron*, 56(1):171–184, 2007.

¹³Z. Füredi and J. Komlos. The eigenvalues of random symmetric matrices. *Combinatorica*, 1(3):233–241, 1981.

¹⁴R.J. Gardner, E. Hermansen, and M. Pachitariu. Toroidal topology of population activity in grid cells. *Nature*, 602:123–128, 2022.

¹⁵W. Gerstner and W. Kistler. *Spiking Neuron Models*. Cambridge University Press, Cambridge, 2002.

¹⁶J. González, M. Cavelli, A. Mondino, N. Rubido, A. Bl Tort, and P. Torterolo. Communication through coherence by means of cross-frequency coupling. *Neuroscience*, 449:157–164, 2020.

¹⁷A. Guilonnet. Bernoulli random matrices. submitted, 2021.

¹⁸A.R. Haig, E. Gordon, J.J. Wright, R.A. Meares, and H. Bahramali. Synchro nous cortical gamma-band activity in task-relevant cognition. *Neuroreport*, 11:669–675, 2000.

¹⁹H. Haken. *Brain Dynamics*. Springer, Berlin, 2002.

²⁰M. Halgren, I. Ulbert, H. Bastuji, D. Fabó, L. Erőss, M. Rey, O. Devinsky, W.K. Doyle, R. Mak-McCully, E. Halgren, L. Wittner, P. Chauvel, G. Heit, E. Eskandar, A. Mandell, and S.S. Cash. The generation and propagation of the human alpha rhythm. *Proc. Natl. Acad. Sci. USA*, 116(47):23772–23782, 2019.

²¹D. Hansel and H. Solominsky. Modeling feature selectivity in local cortical circuits. In C.Koch and I.Segev, editors, *Methods in Neuronal Modeling: From Ions to Networks*, chapter 13. MIT Press, Cambridge, 2 edition, 1998.

²²M. Hashemi, A. Hutt, and J. Sleigh. How the cortico-thalamic feedback affects the EEG power spectrum over frontal and occipital regions during propofol-induced anaesthetic sedation. *J. Comput. Neurosci.*, 39(1):155–159, 2015.

²³P. Holme. Rare and everywhere: Perspectives on scale-free networks. *Nat. Commun.*, 10:1016, 2019.

²⁴X. Huang, W.C. Troy, S.J. Schiff, Q. Yang, H. Ma, C.R. Laing, and J.Y. Wu. Spiral waves in disinhibited mammalian neocortex. *J.Neurosc.*, 24(44):9897–9902, 2004.

²⁵A. Hutt. Brain connectivity reduction reflects disturbed self-organisation of the brain: Neural disorders and general anaesthesia. In V. Cutsuridis, editor, *Multiscale models of Brain Disorders*, pages 207–218. Springer-Verlag, Berlin, 2019.

²⁶A. Hutt. Additive noise-induced system evolution (ANISE). *Front. Appl. Math. Stat.*, 8:879866, 2022.

²⁷A. Hutt and L. Buhry. Study of GABAergic extra-synaptic tonic inhibition in single neurons and neural populations by traversing neural scales: application to propofol-induced anaesthesia. *J. Comput. Neurosci.*, 37(3):417–437, 2014.

²⁸A. Hutt and H. Haken, editors. *Synergetics*. Springer-Verlag, New York, 2020.

²⁹A. Hutt and J. Lefebvre. Additive noise tunes the self-organization in complex systems. In A. Hutt and H. Haken, editors, *Synergetics*, Encyclopedia of Complexity and Systems Science Series, pages 183–196. Springer, New York, 2020.

³⁰A. Hutt and J. Lefebvre. Arousal fluctuations govern oscillatory transitions between dominant γ and α occipital activity during eyes open/closed conditions. *Brain Topography*, 35:108–120, 2021.

³¹A. Hutt, J. Lefebvre, D. Hight, and H. Kaiser. Phase coherence induced by additive gaussian and non-gaussian noise in excitable networks with application to burst suppression-like brain signals. *Front. Appl. Math. Stat.*, 5:69, 2020.

³²A. Hutt, J. Lefebvre, D. Hight, and J. Sleigh. Suppression of underlying neuronal fluctuations mediates EEG slowing during general anaesthesia. *Neuroimage*, 179:414–428, 2018.

³³A. Hutt, A. Mierau, and J. Lefebvre. Dynamic control of synchronous activity in networks of spiking neurons. *PLoS One*, 11(9):e0161488, 2016.

³⁴A. Hutt, S. Rich, T.A. Valiante, and J. Lefebvre. Intrinsic neural diversity quenches the dynamic volatility of neural networks. *Proc. Nat. Acad. Sci. USA*, 120(28):e2218841120, 2023.

³⁵A. Hutt, T. Wahl, N. Voges, J. Hausmann, and J. Lefebvre. Coherence resonance in random erdós-renyi neural networks : mean-field theory. *Front. Appl. Math. Stat.*, 7:697904, 2021.

³⁶M. Iijima, R. Mase, M. Osawa, S. Shimizu, and S. Uchiyama. Event-related synchronization and desynchronization of high-frequency electroencephalographic activity during a visual Go/No-Go paradigm. *Neurophysiobiology*, 71(1):17–24, 2015.

³⁷M. Kaiser, C. C. Hilgetag, and R. Kötter. Hierarchy and dynamics of neural networks. *Front. Neuroinf.*, 4:112, 2010.

³⁸M. Kaschube, M. Schnabel, and F. Wolf. Self-organization and the selection of pinwheel density in visual cortical development. *New J. Phys.*, 10:015009, 2008.

³⁹J.A.S. Kelso. *Dynamic Patterns: The Self-Organization of Brain and Behavior*. MIT Press, Cambridge, 1995.

⁴⁰W. Klimesch, P. Sauseng, and S. Hanslmayr. EEG alpha oscillations: the inhibition-timing hypothesis. *Brain Res. Rev.*, 53(1):63–88, 2007.

⁴¹U. Lee, G. Oh, S. Kim, G. Noh, B. Choi, and G.A. Mashour. Brain networks maintain a scale-free organization across consciousness, anesthesia, and recovery: evidence for adaptive reconfiguration. *Anesthesiology*, 113(5):1081–1091, 2010.

⁴²J. Lefebvre, A. Hutt, and F. Frohlich. Stochastic resonance mediates the state-dependent effect of periodic stimulation on cortical alpha oscillations. *eLife*, 6:e32054, 2017.

⁴³R. M. May. *Stability and complexity in model ecosystems*. Princeton University Press, 2001.

⁴⁴A. Mierau, W. Klimesch, and J. Lefebvre. State-dependent alpha peak frequency shifts: Experimental evidence, potential mechanisms and functional implications. *Neuroscience*, 360:146–154, 2017.

⁴⁵V.B. Mountcastle. Modality and topographic properties of single neurons of cat's somatic sensory cortex. *Neurophysiol.*, 20:408

⁴⁶G. Nicolis and I. Prigogine. *Self-Organization in Non-Equilibrium Systems: From Dissipative Structures to Order Through Fluctuations*. J. Wiley and Sons, New York, 1977.

⁴⁷G. Pfurtscheller and F.H. Lopes da Silva. Event-related EEG/MEG synchronization and desynchronization: basic principles. *Clin. Neurophysiol.*, 110(11):1842–1857, 1999.

⁴⁸B. Pietras, V. Schmutz, and T. Schwalger. Mesoscopic description of hippocampal replay and metastability in spiking neural networks with short-term plasticity. *PLoS Comput Biol*, 18(12):e1010809, 2022.

⁴⁹A.S. Pikovsky and J. Kurths. Coherence resonance in a noise-driven excitable system. *Phys. Rev. Lett.*, 78:775–778, 1997.
<footnote n="50">⁵⁰J.C. Prechtl, L.B. Cohen, B. Pesaran, P.P. Mitra, and D. Kleinfeld. Visual stimuli induce waves of electrical activity in turtle cortex. *Proc. Nat. Acad. Sci.*, 94:7621–7626, 1997.</footnote>

<footnote n="51">⁵¹S. Ray and J.H. Maunsell. Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. *PLoS Biol.*, 9(4):e1000610, 2011.</footnote>

<footnote n="52">⁵²S. Rich, H. M. Chameh, J. Lefebvre, and T. A. Valiante. Resilience through diversity: Loss of neuronal heterogeneity in epileptogenic human tissue impairs network resilience to sudden changes in synchrony. *Cell Rep.*, 39(8):110863, 2022.</footnote>

<footnote n="53">⁵³S. Rich, A. Hutt, F.K. Skinner, T.A. Valiante, and J. Lefebvre. Neurostimulation stabilizes spiking neural networks by disrupting seizure-like oscillatory transitions. *Sci. Rep.*, 10:15408, 2020.</footnote>

<footnote n="54">⁵⁴J. Riedinger and A. Hutt. Mathematical model insights into eeg origin under transcranial direct current stimulation (tdcs) in the context of psychosis. *J. Clin. Med.*, 11(7):1845, 2022.</footnote>

<footnote n="55">⁵⁵P.A. Robinson, P.N. Loxley, S.C. O'Connor, and C.J. Rennie. Modal analysis of corticothalamic dynamics, electroencephalographic spectra and evoked potentials. *Phys. Rev. E*, 63:041909, 2001.</footnote>

<footnote n="56">⁵⁶M. Seeber, R. Scherer, J. Wagner, T. Solis-Escalante, and G.R. Müller-Putz. Eeg beta suppression and low gamma modulation are different elements of human upright walking. *Front. Hum. Neurosci.*, 8:485, 2014.</footnote>

<footnote n="57">⁵⁷T. Shibata, I. Shimoyama, T. Ito, D. Abla, H. Iwasa, K. Koseki, N. Yamanouchi, T. Sato, and Y. Nakajima. Event-related dynamics of the gamma-band oscillation in the human brain: information processing during</footnote>

<footnote>a Go/No-Go hand movement task. *Neurosci. Res.*, 33(3):215–222, 1999.</footnote>

<footnote n="58">⁵⁸L. Solanka, M.C.W. van Rossum, and M.F. Nolan. Noise promotes independent control of gamma oscillations and grid firing within recurrent attractor networks. *eLife*, 4:e06444, 2015.</footnote>

<footnote n="59">⁵⁹M. Takemi, Y. Masakado, M. Liu, and J. Ushiba. Event-related desynchronization reflects downregulation of intracortical inhibition in human primary motor cortex. *Journal of Neurophysiology*, 110(5):1158–1166, 2013.</footnote>

<footnote n="60">⁶⁰S. Vijayan and N.J. Kopell. Thalamic model of awake alpha oscillations and implications for stimulus processing. *Proc. Natl. Acad. Sci. USA*, 109(45):18553–18558, 2012.</footnote>

<footnote n="61">⁶¹M. A. Whittington, M. O. Cunningham, F.E.N. Le Beau, C. Racca, and R.D. Traub. Multiple origins of the cortical gamma rhythm. *Developmental Neurobiology*, 71:92–106, 2011.</footnote>

<footnote n="62">⁶²E.P. Wigner. On the distribution of the roots of certain symmetric matrices. *Ann. Math.*, 67:325–327, 1958.</footnote>

<footnote n="63">⁶³T. W. Wilson, T. J. McDermott, M. S. Mills, N. M. Coolidge, and E. Heinrichs-Graham. tDCS modulates visual gamma oscillations and basal alpha activity in occipital cortices: Evidence from MEG. *Cerebral Cortex*, 28:1597–1609, 2018.</footnote>

<footnote n="64">⁶⁴P. Zhou, S.D. Burton, N.N. Urban, and G.B. Ermentrout. Impact of neuronal heterogeneity on correlated colored noise-induced synchronization. *Front Comput Neurosci*, 7:113, 2013.</footnote>
<footnote n="50">⁵⁰J.C. Prechtl, L.B. Cohen, B. Pesaran, P.P. Mitra, and D. Kleinfeld. Visual stimuli induce waves of electrical activity in turtle cortex. *Proc. Nat. Acad. Sci.*, 94:7621–7626, 1997.</footnote>

<footnote n="51">⁵¹S. Ray and J.H. Maunsell. Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. *PLoS Biol.*, 9(4):e1000610, 2011.</footnote>

<footnote n="52">⁵²S. Rich, H. M. Chameh, J. Lefebvre, and T. A. Valiante. Resilience through diversity: Loss of neuronal heterogeneity in epileptogenic human tissue impairs network resilience to sudden changes in synchrony. *Cell Rep.*, 39(8):110863, 2022.</footnote>

<footnote n="53">⁵³S. Rich, A. Hutt, F.K. Skinner, T.A. Valiante, and J. Lefebvre. Neurostimulation stabilizes spiking neural networks by disrupting seizure-like oscillatory transitions. *Sci. Rep.*, 10:15408, 2020.</footnote>

<footnote n="54">⁵⁴J. Riedinger and A. Hutt. Mathematical model insights into eeg origin under transcranial direct current stimulation (tdcs) in the context of psychosis. *J. Clin. Med.*, 11(7):1845, 2022.</footnote>

<footnote n="55">⁵⁵P.A. Robinson, P.N. Loxley, S.C. O'Connor, and C.J. Rennie. Modal analysis of corticothalamic dynamics, electroencephalographic spectra and evoked potentials. *Phys. Rev. E*, 63:041909, 2001.</footnote>

<footnote n="56">⁵⁶M. Seeber, R. Scherer, J. Wagner, T. Solis-Escalante, and G.R. Müller-Putz. Eeg beta suppression and low gamma modulation are different elements of human upright walking. *Front. Hum. Neurosci.*, 8:485, 2014.</footnote>

<footnote n="57">⁵⁷T. Shibata, I. Shimoyama, T. Ito, D. Abla, H. Iwasa, K. Koseki, N. Yamanouchi, T. Sato, and Y. Nakajima. Event-related dynamics of the gamma-band oscillation in the human brain: information processing during</footnote>

<footnote>a Go/No-Go hand movement task. *Neurosci. Res.*, 33(3):215–222, 1999.</footnote>

<footnote n="58">⁵⁸L. Solanka, M.C.W. van Rossum, and M.F. Nolan. Noise promotes independent control of gamma oscillations and grid firing within recurrent attractor networks. *eLife*, 4:e06444, 2015.</footnote>

<footnote n="59">⁵⁹M. Takemi, Y. Masakado, M. Liu, and J. Ushiba. Event-related desynchronization reflects downregulation of intracortical inhibition in human primary motor cortex. *Journal of Neurophysiology*, 110(5):1158–1166, 2013.</footnote>

<footnote n="60">⁶⁰S. Vijayan and N.J. Kopell. Thalamic model of awake alpha oscillations and implications for stimulus processing. *Proc. Natl. Acad. Sci. USA*, 109(45):18553–18558, 2012.</footnote>

<footnote n="61">⁶¹M. A. Whittington, M. O. Cunningham, F.E.N. Le Beau, C. Racca, and R.D. Traub. Multiple origins of the cortical gamma rhythm. *Developmental Neurobiology*, 71:92–106, 2011.</footnote>

<footnote n="62">⁶²E.P. Wigner. On the distribution of the roots of certain symmetric matrices. *Ann. Math.*, 67:325–327, 1958.</footnote>

<footnote n="63">⁶³T. W. Wilson, T. J. McDermott, M. S. Mills, N. M. Coolidge, and E. Heinrichs-Graham. tDCS modulates visual gamma oscillations and basal alpha activity in occipital cortices: Evidence from MEG. *Cerebral Cortex*, 28:1597–1609, 2018.</footnote>

<footnote n="64">⁶⁴P. Zhou, S.D. Burton, N.N. Urban, and G.B. Ermentrout. Impact of neuronal heterogeneity on correlated colored noise-induced synchronization. *Front Comput Neurosci*, 7:113, 2013.</footnote>