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

HAL
open science

# Discrete Synaptic Events Induce Global Oscillations in Balanced Neural Networks

Denis S Goldobin, Matteo Di Volo, Alessandro Torcini

► To cite this version:

Denis S Goldobin, Matteo Di Volo, Alessandro Torcini. Discrete Synaptic Events Induce Global Oscillations in Balanced Neural Networks. Physical Review Letters, 2024, 133 (23), pp.238401. 10.1103/PhysRevLett.133.238401. hal-04943715

HAL Id: hal-04943715

https://hal.science/hal-04943715v1

Submitted on 12 Feb 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Discrete synaptic events induce global oscillations in balanced neural networks

Denis S. Goldobin,¹,² Matteo di Volo,³ and Alessandro Torcini⁴,⁵,⁶,*

¹Institute of Continuous Media Mechanics, Ural Branch of RAS, Acad. Korolev street 1, 614013 Perm, Russia

²Institute of Physics and Mathematics, Perm State University, Bukirev street 15, 614990 Perm, Russia

³Université Claude Bernard Lyon 1, Institut National de la Santé et de la Recherche Médicale,
Stem Cell and Brain Research Institute U1208, Bron, France

⁴Laboratoire de Physique Théorique et Modélisation, CY Cergy Paris Université,
CNRS, UMR 8089, 95302 Cergy-Pontoise cedex, France

⁵CNR - Consiglio Nazionale delle Ricerche - Istituto dei Sistemi Complessi,
via Madonna del Piano 10, 50019 Sesto Fiorentino, Italy

⁶INFN Sezione di Firenze, Via Sansone 1, 50019 Sesto Fiorentino, Italy

(Dated: November 1, 2024)

Despite neural dynamics is triggered by discrete synaptic events, the neural response is usually obtained within the diffusion approximation (DA) representing the synaptic inputs as Gaussian noise. We derive a mean-field formalism encompassing synaptic shot-noise for sparse balanced neural networks. For low (high) excitatory drive (inhibitory feedback) global oscillations (GOs) emerge via continuous or hysteretic transitions, correctly predicted by our approach, but not from the DA. At sufficiently low in-degrees the nature of these GOs changes from drift-driven to cluster activation.

*Introduction.* In several contexts the discrete nature of stochastic events should be taken into account to correctly predict the system dynamics. A typical example is represented by shot-noise, which is conveyed by pulses and is therefore discontinous, at variance with white noise, which is associated to thermal fluctuations and is continuous [1]. The inclusion of shot-noise is fundamental to fully characterize the emergent phenomena in many fields of physics ranging from mesoscopic conductors [2] to driven granular gases [3].

The discrete nature of the events is an innate characteristic also of the neural dynamics, where a neuron receives inputs from other neurons via electrical pulses, that manifest as post-synaptic potentials (PSPs). The PSPs stimulating a neuron in the cortex are usually assumed to be uncorrelated with small amplitudes and high arrival rates. Therefore the mean-field (MF) neural dynamics has been examined within the framework of the Diffusion Approximation (DA) [4, 5] by treating synaptic inputs as a continuous Gaussian process.

However, several experiments have shown that rare PSPs of large amplitude can have a fundamental impact on the cortical activity [6, 7] and that synaptic weight distributions display a long tail towards large amplitudes [8–10]. Furthermore, networks of inhibitory neurons with low connectivity (in-degree $K \approx 30–80$) have been identified in the cat visual cortex [11] and in the rat hippocampus [12], where they are believed to be at the origin of global oscillations (GOs) in the $\gamma$-band [13]. Moreover, the cortical connectivity is definitely more sparse in primates than in rodents as recently shown [14].

These experimental evidences call for the development of a MF formalism able to incorporate the effect of discrete synaptic events for diluted random neural networks. Population based formalisms, taking into account synaptic shot-noise, have been previously developed for

Integrate-and-Fire models [15–18]. However, such approaches are limited to stationary solutions and cannot describe the emergence of oscillatory behaviours.

In this Letter, we introduce a complete mean-field (CMF) approach for balanced neural networks [19, 20], taking into account the sparseness of the network and the discreteness of the synaptic pulses, able to reproduce all the possible dynamical states. For simplicity, but without any loss of generality, we consider inhibitory balanced networks subject to an external excitatory drive [21–24].

Firstly, we illustrate that the DA fails in reproducing oscillatory dynamics in spiking neural networks for sufficiently low excitatory drive (high inhibitory feedback) by considering conductance- and current-based models. However, this regime is correctly reproduced by a MF approach whenever the sparse and discrete synaptic inputs are taken into account. Moreover, for Quadratic Integrate-and-Fire (QIF) [25, 26] networks via the CMF approach we obtain a complete bifurcation diagram encompassing asynchronous (ARs) and oscillatory regimes (ORs), where individual neurons spike irregularly. The CMF reveals sub- and super-critical Hopf bifurcations from the AR to the OR as well as a region of coexistence of these two phases not captured by the DA [27]. Event-driven simulations of large QIF networks confirm the scenario predicted within the CMF theory. Furthermore, we show that the GOs, induced by discrete synaptic events, can emerge due to two different mechanisms: cluster activation at sufficiently small $K$ and drift-driven at larger $K$.

*The balanced network.* As a prototype of a dynamically balanced system we consider a sparse inhibitory network made of $N$ pulse-coupled neurons whose mem-
brane potentials evolve according to the equations

$$ \dot{v}_i(t) = F(v_i) + I - g \sum_{j=1}^{N} \sum_{n} \epsilon_{ji} \delta(t - t_j^{(n)}) ; \quad (1) $$

where $I$ represents an external DC current, $g$ the synaptic coupling, and the last term the inhibitory synaptic current. The latter is the linear superposition of instantaneous inhibitory PSPs emitted at times $t_j^{(n)}$ from the pre-synaptic neurons connected to neuron $i$. $\epsilon_{ji}$ is the adjacency matrix of the random network with entries 1 (0) if the connection from node $j$ to $i$ exists (or not), and we assume the same in-degree $K = \sum_j \epsilon_{ji}$ for all neurons. We consider two paradigmatic models of spiking neuron: the quadratic integrate-and-fire (QIF) with $F(v) = v^2$ [23–25, 28, 29], which is a current-based model of class I excitability; and the Morris-Lecar (ML) [30], a conductance-based model representing a class II excitable membrane [31]. The DC current and the synaptic coupling are assumed to scale as $I = i_0 \sqrt{K}$ and $g = g_0 / \sqrt{K}$ as usually done in order to ensure a self-sustained balanced state for sufficiently large $K$ [19, 22–24, 32, 33].

**Mean-field description.** For a sufficiently sparse network with $K \ll N$, the spike trains emitted by $K$ presynaptic neurons can be assumed to be uncorrelated and Poissonian [20, 21], therefore the MF dynamics of a generic neuron can be represented in terms of the following Langevin equation:

$$ \dot{V}(t) = F(V) + I - gS(t) \quad (2) $$

where $S(t)$ is a Poissonian train of $\delta$-spikes with rate $R(t) = K\nu(t)$, and $\nu(t)$ the population firing rate self-consistently estimated. Usually the Poissonian spike trains are approximated within the DA [5, 34] as $S(t) = R(t) + \sqrt{R(t)}\xi(t)$, where $\xi(t)$ is a Gaussian white noise term. However, the DA can fail in reproducing the neural dynamics. Indeed, as shown in Fig. 1 (a) for a sparse ML network, by employing the DA in (2) one obtains an asynchronous dynamics (blue curve), while the network evolution, characterized by GOs with frequency $f_C \approx 18$ Hz (black dots), can be recovered only by explicitly taking into account the Poissonian spike trains in (2) (red line).

In the MF framework the population dynamics is usually described in terms of the membrane potential probability distribution function $P(V, t)$, whose time evolution for the QIF model is given (according to (2)) by the continuity equation

$$ \partial_t P(V, t) + \partial_V [(V^2 + I)P(V, t)] = R(t)\Delta P(V, T) \quad (3) $$

with boundary condition $\lim_{V \to \infty} V^2 P(V, t) = \nu(t)$ and where $\Delta P(V, T) = [P(V^+, t) - P(V, t)]$ with $V^+ = V + g$. By assuming that $g$ is sufficiently small we can expand the latter term as $\Delta P(V, t) = \sum_{p=1}^{\infty} \frac{g_p}{p!} \partial_V^p P(V, t)$; and by limiting to the first two terms in this expansion we

<figure><img src="image_3.png" /><figcaption>FIG. 1. Population firing rate $\nu(t)$ versus time for ML (a) and QIF (b) models. Black circles refer to network simulations, blue (red) solid line to diffusive (shot-noise) Langevin results obtained by integrating Eq. (2). In (b) magenta (green) dashed line denotes the DA (shot-noise) MF approximation resulting from the integration of Eq. (3) (Eq. (4)) for the QIF model, more details in [31]. The parameters for the ML model are $K = 20$, $i_0 = 0.1$, $g_0 = 5$ and network size $N = 20000$, other parameters are reported in the supplemental material [31]. For the QIF model $K = 200$, $i_0 = 0.16$, $g_0 = 4$ and $N = 80000$ [36].</figcaption></figure>

recover the DA corresponding to the following Fokker-Planck Equation (FPE) [35]

$$ \partial_t P(V, t) + \partial_V [(V^2 + A(t))P(V, t)] = D(t)\partial_V^2 P(V, t) \quad (4) $$

where $A(t) = \sqrt{K}[i_0 - g_0\nu(t)]$ and $D(t) = g_0^2\nu(t)/2$. The DA can give incorrect predictions for the QIF model, as well. Indeed as shown in Fig. 1 (b) the network dynamics is oscillatory with $f_c \approx 40$ Hz (black circles). This evolution is correctly captured by the MF equation (3) (green dashed line) and by the Langevin equation (2) driven by shot-noise (red solid line). On the contrary the FPE (4) (dashed magenta line) and the diffusive Langevin formulation (blue solid line) converge to a stable fixed point corresponding to asynchronous dynamics. Therefore to reproduce the collective dynamical regimes observable in the network it is necessary to consider the continuity equation (3). In this respect we have developed a CMF formalism encompassing synaptic shot-noise to identify the various possible regimes displayed by (3) and to analyse their stability.

The QIF model evolution can be transformed in that of a phase oscillator, the so-called $\theta$-neuron [25, 37], by introducing the phase variable $\theta = 2 \arctan V$. However, this transformation has the drawback that even uncoupled neurons are associated to a non flat distribution of the phases, thus rendering quite difficult to distinguish asynchronous from partially synchronized regimes [38, 39]. A more appropriate phase transformation to analyse the synchronization phenomena is the following $\psi = 2 \arctan (V/\sqrt{I}) \in [-\pi, \pi]$, which leads to a uniformly rotating phase in the absence of incoming pulses for supra-threshold neurons with $I > 0$ [31].

By considering the probability distribution of the phases $w(\psi, t) = P(V, t)(I + V^2)/(2\sqrt{I})$, Eq. (3) can be rewritten in terms of the so-called Kuramoto–Daido order parameters $z_n$ [40, 41] by expanding in Fourier space the
distribution as $w(\psi, t) = (2\pi)^{-1} \sum_{n=-\infty}^{+\infty} z_n e^{-in\psi}$ with $z_0 = 1$ and $z_{-n} = z_n^*$. After laborious but straightforward calculations, one obtains the following evolution equations

$$ \dot{z}_n = i2n\sqrt{I}z_n + K\nu \left[ \sum_{m=0}^{+\infty} I_{nm}(\alpha) z_m - z_n \right], \quad (5) $$

where $n = 1, 2, 3, ...$, $\alpha \equiv g/\sqrt{I} = g_0/(\sqrt{i_0}K^{3/4})$, and the explicit expressions for $I_{nm}(\alpha)$ are reported in [31].

The firing rate can be self-consistently determined by the flux at the firing threshold $\lim_{V \to \infty} V^2 P(V, t) = 2\sqrt{I}w(\pi, t)$, as follows

$$ \nu = 2\sqrt{I}w(\pi, t) = \frac{\sqrt{I}}{\pi} \text{Re} \left( 1 + 2 \sum_{k=1}^{\infty} (-1)^k z_k \right). \quad (6) $$

The dynamics of the system (5,6) is controlled by only two parameters: $K$ and $\alpha$. Thus, we can limit to derive a bidimensional phase diagram in the plane $(i_0/g_0^2, K)$, that will comprehensively cover the entire diversity of the macroscopic regimes observable in the network. In particular, we have estimated the stationary solutions of Eqs. (5,6) by truncating the Fourier expansion in (5) to $M \ge 100$ modes in order to guarantee a numerical accuracy of $O(10^{-12})$ for all the parameter values. The linear stability of the asynchronous state joined to the derivation of the corresponding amplitude equations (via a weakly nonlinear approach) has allowed us to identify the HB line where the oscillatory dynamics emerges together with the super- or sub-critical nature of the bifurcations (for more details see section S3 in [31]). The HB line obtained via the CMF (within the DA) is reported as an orange (black) line in Fig. 2 (a). While the HBs are always super-critical within the DA, the HBs induced by the shot-noise can be either super- (solid orange line) or sub-critical (dashed orange line), thus allowing for regions where AR and OR coexist, see Fig. 2 (b). A peculiarity of the CMF results is that the HB line is reentrant, thus in a certain range of $i_0/g_0^2$ we have an AR only within a finite interval of in-degrees and GOs at sufficiently small and large $K$ (as shown in Fig. 2 (c)). As explained in the following these two oscillatory regimes are due to different mechanisms.

Furthermore, there is a dramatic difference among the 2 MF approaches at small $i_0$ (large $g_0$): within the DA GOs are observable only above a critical $K$ diverging to infinity for $i_0/g_0^2 \to 0$, while for the CMF analysis GOs are present at any $K$ value for $i_0/g_0^2 < 0.00029$.

**Network Simulations.** In order to verify the CMF predictions we have performed essentially exact numerical simulations of QIF networks, according to (1), by employing an event-driven integration scheme [42], which allowed us to follow the network dynamics for long times, up to 50–100 sec, for system of sizes $N = 10000-80000$

<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Phase diagram for the QIF network in the plane $(i_0/g_0^2, K)$: the black solid line is the super-critical HB line obtained within the DA; the orange solid (dashed) line is the super- (sub-) critical HB line given by the CMF; the symbols refer to numerical estimations of the HBs and SNBs. The green (blue) circles denote HBs obtained by performing quasi-adiabatic simulations by varying $K$ ($i_0$) for constant $i_0$ ($K$) values; the magenta stars indicate SNBs. For more details see [31]. (b-c) Average order parameter $\rho_A$ versus $i_0$ ($K$) for quasi-adiabatic simulations: black circles refer to decreasing (increasing) $i_0$ ($K$), while red ones to increasing $i_0$. The blue dashed line in (b) denotes the sub-critical HB given by the CMF and the magenta dot-dashed line to numerically estimated SNB; the two green dashed lines in (c) indicate the HBs given by the CMF. In the insets in (b-c) the population firing rates $\nu(t)$ versus time are reported for the states indicated by the corresponding colored filled circles, in the insets in (b) the results for the CMF (5) are also shown as green solid lines. The values of $\rho_A$ in (b) ((c)) refer to $K = 100$ ($i_0 = 0.00055$) averaged over 5 network realizations, with $N = 80000$, for 30 s following a transient of 20 s. For all data $g_0 = 1$ [36].</figcaption></figure>

[31]. In particular, to characterize the macroscopic evolution of the network we measured the indicator introduced in [43] $\rho = (\sigma_V^2/\sigma_i^2)^{1/2}$, where $\sigma_i^2 = \langle v_i^2 \rangle - \langle v_i \rangle^2$ and $\sigma_V$ is the standard deviation of the mean membrane potential $\bar{V} = \sum_{i=1}^N v_i/N$, with $\cdot$ denoting an ensemble (a time) average. A coherent (asynchronous) macroscopic activity is characterized by a value of $\rho$ remaining finite (vanishing as $\rho \propto N^{-1/2}$) for $N \to \infty$ [24, 44]. The actual value of $\rho$ is related to the level of synchronization among the neurons: perfect synchrony to $\rho \equiv 1$.

A finite size analysis of the order parameter $\rho_A$ averaged over several different network realizations has allowed us to identify the HBs and the Saddle-Node Bifurcations (SNBs) of limit cycles displayed in Fig. 2 [31]. In particular, in Fig. 2 (a) green (blue) circles refer to HBs identified via quasi-adiabatic simulations by varying $K$ ($i_0$) for constant $i_0$ ($K$) values; while the magenta stars indicate SNBs. Numerical simulations are in good agreement with the CMF results and allowed us also the iden-
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a-b) (Lower Panels) Effective input current $A(t)$ versus time: red lines (black dots) refer to CMF (network simulations) results. (Upper Panels) Membrane potential evolution in time: black lines (other colors) refer to $V_0(t)$ (single neuron dynamics). (c) Frequency of the GOs $f_c$ (black circles) versus $K$, the red solid line refers to $f_0$, green vertical dashed lines have the same meaning as in Fig. 2 (c). Data correspond to $i_0/g_0^2 = 0.00055$ (with $g_0 = 20$), to $K = 10$ (a) and $K = 210$ (b), network simulations to $N = 80000$ [36].</figcaption></figure>

tification of a coexistence region for asynchronous and oscillatory collective dynamics. A hysteretic transition from AR to OR obtained by varying quasi-adiabatically $i_0$ is displayed in Fig. 2 (b), the coexistence region can be clearly identified between the sub-critical HB (blue dashed line) and the SNB (magenta dashed line). Two coexisting solutions are reported in the insets of Fig. 2 (b) confirming the good agreement between CMF (green lines) and the network simulations (red and black lines). Furthermore, as shown in Fig. 2 (c) for sufficiently small $i_0/g_0^2$ values GOs are observable at small ($K \le 30$) and large ($K \ge 200$) in-degrees, while the AR is present only at intermediate in-degrees ($K \in [40:180]$). The dynamics in these 3 intervals is visualized by reporting in the insets of Fig. 2 (c) the firing rates $\nu(t)$ at $K = 10$ (black line), $K = 60$ (blue line) and $K = 210$ (green line).

In large part of the phase diagram (namely, for $i_0/g_0^2 < 0.2$), both in the AR and OR we observe an irregular firing activity of the neurons associated to mean coefficient of variations $CV \simeq O(1)$ [45], as expected in sparse balanced networks.

**Two kinds of GOs.** As previously mentioned, we can identify two classes of GOs induced by discrete synaptic events in the interval $i_0/g_0^2 \in [0.00036:0.00070]$. Their difference is already clear by considering the MF membrane potential evolution $V_0(t)$ given by the following zeroth-order Langevin equation for the QIF:

$$ \dot{V}_0(t) = V_0^2 + \sqrt{K}[i_0 - g_0\nu(t)] = V_0^2 + A(t) \quad ; \quad (7) $$

where current fluctuations have been neglected and $\nu(t)$ is the population firing induced by the shot-noise. Whenever $A < 0$ ($A > 0$) the QIF model displays excitable dynamics (periodic firing) [26]. The GOs reported in the insets of Fig. 2 (c) for $K = 10$ ($K = 210$) are characterized

by $A(t)$ always negative (positive for large part of the oscillation period) as shown in Fig. 3 (a) (Fig. 3 (b)) (lower panels). Therefore, for $K = 10$ ($K = 210$) $V_0(t)$ displays sub-threshold oscillations (large excursions from negative to positive values driven by $A(t) > 0$) as shown in the upper panel of Fig. 3 (a) (Fig. 3 (b)).

As shown in Fig. 2 (c), the low (high) in-degree GOs are characterized by a low (high) level of coherence among the neurons, this is confirmed by the evolution of the membrane potentials $v_i(t)$ of 4 generic neurons reported in Fig. 3 (a-b). In both cases the neurons spike irregularly [46], however for $K = 210$ the single neurons $v_i(t)$ essentially follow the MF evolution $V_0(t)$, while for $K = 10$ their dynamics is quite uncorrelated.

These behaviours can be explained by 2 different mechanisms once noticed that for both cases the GO frequency $f_C$ is extremely close to the firing frequency of an isolated neuron $f_0 = 1/T_0 = \sqrt{I}/\pi$, see Fig. 3 (c). This suggests that at a first approximation the GOs are due to the neurons not receiving any inhibitory PSP from reset to threshold. For low $K$, whenever a neuron fires large amplitude inhibitory PSPs are delivered. These induce a transient synchronization in the $K$ post-synaptic neurons and a sub-group, not receiving further PSPs, can eventually reach threshold together at a time $\simeq T_0$. This transient synchronizing effect of small clusters of neurons (termed cluster activation [47]) is at the basis of the GOs observable for $K = 10$ in Fig. 2 (c). For increasing $K$, the amplitude of the PSPs decreases, therefore above some critical in-degree ($K \simeq 30$ in this case) a single inhibitory PSP is no more able to induce a sufficiently strong synchronizing effect on the post-synaptic neurons and the dynamics becomes asynchronous (as shown in Fig. 2 (c)).

For larger $K$, the post-synaptic neurons receive several small inhibitory PSPs at each population burst, whenever $K$ is sufficiently large a non negligible part of the neurons can get synchronized by the discharge of inhibitory PSPs. As shown in Fig 3 (b), the time courses of the membrane potentials are now extremely coherent by approaching the threshold, where fluctuations lead to irregular firing of the neurons. However, a sufficient percentage of neurons *drift-driven* is always able to fire together with a period $\simeq T_0$ giving rise to the GOs.

**Conclusions.** We have shown that the macroscopic phase-diagram of balanced networks is strongly influenced by the discreteness and the finite amplitude of PSPs. In particular, we have developed a CMF formalism by including Poissonian shot-noise which reproduces quite well the network simulations, at variance with the DA. This scenario is robust and extends beyond instantaneous synapses to exponentially decaying PSPs, as shown in Fig. S4 in [31].

Our analysis of balanced inhibitory networks has revealed the existence of two kinds of GOs induced by discrete synaptic events, thus completing the previous scenario
nario based on the DA [21, 27]. Furthermore, we have shown that GOs can emerge even in extremely sparse inhibitory networks with frequencies going from 1-2 Hz (δ-band) to 100 Hz (γ-band). Thus providing theoretical support for the supposition reported in [13] that the γ-oscillations observed in the hippocampus are generated by sub-networks of interneurons with low in-degrees $K \approx 30-80$ [12].

The CMF approach is valid in sparse networks for $K \ll N$. Whenever $K \approx O(N)$ the correlations among the spike trains reaching the neuron should be taken in account, an extension of the CMF in this direction will be worth of future investigations. The effect of finite $N$ fluctuations has been analyzed in globally coupled QIF networks [48, 49], it will be interesting to extend such approach to random networks. Finally, the CMF formalism can be generalised to neural systems with delay and synaptic kinetics as shown in sub-section S3.G in [31], this will be the subject of future studies.

We acknowledge stimulating discussions with Alberto Bacci, Alberto Ferrara, Nina La Miciotta, Lyudmila Klimenko, Gianluigi Mongillo, Simona Olmi, Antonio Politi. We are also grateful to the anonymous Referees for extremely useful suggestions. D.S.G. acknowledges the support of the CNR Short Term Mobility Programme 2021 for a visit to Istituto dei Sistemi Complessi, Sesto Fiorentino, Italy where part of this work was developed. A.T. received financial support by the Labex MME-DII (Grant No. ANR-11-LBX-0023-01), by CY Generations (Grant No ANR-21-EXES-0008), and together with M.V. by the ANR Project ERMUNDY (Grant No. ANR-18-CE37-0014) and M.V. by the Labex CORTEX (Grant No. ANR-11-LABX-0042) of Université Claude Bernard Lyon 1 and by the the ANR via the Junior Professor Chair in Computational Neurosciences Lyon 1.

<footnote>* corresponding author: alessandro.torcini@cyu.fr</footnote>

<footnote>[1] W. Schottky, Annalen der Physik **362**, 541 (1918).</footnote>

<footnote>[2] Y. M. Blanter and M. Büttiker, Physics Reports **336**, 1 (2000).</footnote>

<footnote>[3] D. Lucente, M. Viale, A. Gnoli, A. Puglisi, and A. Vulpiani, Physical Review Letters **131**, 078201 (2023).</footnote>

<footnote>[4] R. Capocelli and L. Ricciardi, Kybernetik **8**, 214 (1971).</footnote>

<footnote>[5] H. C. Tuckwell, *Introduction to theoretical neurobiology: nonlinear and stochastic theories*, Vol. 2 (Cambridge University Press, 1988).</footnote>

<footnote>[6] S. Song, P. J. Sjöström, M. Reigl, S. Nelson, and D. B. Chklovskii, PLoS Biology **3**, e68 (2005).</footnote>

<footnote>[7] S. Lefort, C. Tomm, J.-C. F. Sarria, and C. C. Petersen, Neuron **61**, 301 (2009).</footnote>

<footnote>[8] R. Miles, The Journal of Physiology **431**, 659 (1990).</footnote>

<footnote>[9] B. Barbour, N. Brunel, V. Hakim, and J.-P. Nadal, TRENDS in Neurosciences **30**, 622 (2007).</footnote>

<footnote>[10] G. Buzsáki and K. Mizuseki, Nature Reviews Neuroscience **15**, 264 (2014).</footnote>

<footnote>[11] Z. F. Kisvárday, C. Beaulieu, and U. T. Eysel, Journal of Comparative Neurology **327**, 398 (1993).</footnote>

<footnote>[12] A. Sik, M. Penttonen, A. Ylinen, and G. Buzsáki, Journal of Neuroscience **15**, 6651 (1995).</footnote>

<footnote>[13] G. Buzsáki and X.-J. Wang, Annual review of neuroscience **35**, 203 (2012).</footnote>

<footnote>[14] G. A. Wildenberg, M. R. Rosen, J. Lundell, D. Paukner, D. J. Freedman, and N. Kasthuri, Cell Reports **36** (2021).</footnote>

<footnote>[15] M. J. E. Richardson and R. Swarbrick, Physical Review Letters **105**, 178102 (2010).</footnote>

<footnote>[16] R. Iyer, V. Menon, M. Buice, C. Koch, and S. Mihalas, PLoS Computational Biology **9**, e1003248 (2013).</footnote>

<footnote>[17] S. Olmi, D. Angulo-Garcia, A. Imparato, and A. Torcini, Scientific Reports **7**, 1577 (2017).</footnote>

<footnote>[18] F. Droste and B. Lindner, Journal of Computational Neuroscience **43**, 81 (2017).</footnote>

<footnote>[19] C. van Vreeswijk and H. Sompolinsky, Science **274**, 1724 (1996).</footnote>

<footnote>[20] N. Brunel, Journal of Computational Neuroscience **8**, 183 (2000).</footnote>

<footnote>[21] N. Brunel and V. Hakim, Neural computation **11**, 1621 (1999).</footnote>

<footnote>[22] J. Kadmon and H. Sompolinsky, Phys. Rev. X **5**, 041030 (2015).</footnote>

<footnote>[23] M. Monteforte and F. Wolf, Phys. Rev. Lett. **105**, 268104 (2010).</footnote>

<footnote>[24] M. di Volo and A. Torcini, Phys. Rev. Lett. **121**, 128301 (2018).</footnote>

<footnote>[25] G. B. Ermentrout and N. Kopell, SIAM Journal on Applied Mathematics **46**, 233 (1986).</footnote>

<footnote>[26] B. Gutkin, in *Encyclopedia of computational neuroscience* (Springer, 2022) pp. 3412–3419.</footnote>

<footnote>[27] M. Di Volo, M. Segneri, D. S. Goldobin, A. Politi, and A. Torcini, Chaos: An Interdisciplinary Journal of Nonlinear Science **32**, 023120 (2022).</footnote>

<footnote>[28] C. R. Laing, The Journal of Mathematical Neuroscience **8**, 1 (2018).</footnote>

<footnote>[29] I. Ratas and K. Pyragas, Physical Review E **100**, 052211 (2019).</footnote>

<footnote>[30] C. Morris and H. Lecar, Biophysical journal **35**, 193 (1981).</footnote>

<footnote>[31] See Supplemental Material at [URL will be inserted by publisher], which includes Refs. [24, 25, 37, 40–42, 45, 50, 51], for details on the employed models, on the integration of the neural networks as well as of the population models, and on the CMF approach.</footnote>

<footnote>[32] A. Renart, J. de la Rocha, P. Bartho, L. Hollender, N. Parga, A. Reyes, and K. D. Harris, Science **327**, 587 (2010).</footnote>

<footnote>[33] A. Litwin-Kumar and B. Doiron, Nat Neurosci **15**, 1498 (2012).</footnote>

<footnote>[34] R. Capocelli and L. Ricciardi, Kybernetik **8**, 214 (1971).</footnote>

<footnote>[35] E. Haskell, D. Q. Nykamp, and D. Tranchina, Network: Computation in Neural Systems **12**, 141 (2001).</footnote>

<footnote>[36] The times (frequencies) are reported in physical units by assuming a membrane time constant $\tau_m = 10$ ms.</footnote>

<footnote>[37] B. Ermentrout, Scholarpedia **3**, 1398 (2008), revision #122134.</footnote>

<footnote>[38] B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, Physical Review E **76**, 055201(R) (2007).</footnote>

<footnote>[39] A. V. Dolmatova, D. S. Goldobin, and A. Pikovsky, Physical Review E **96**, 062204 (2017).</footnote>
[40] Y. Kuramoto, *Chemical oscillations, waves, and turbulence*, Vol. 19 (Springer Science & Business Media, 2012).

[41] H. Daido, *Progress of Theoretical Physics* **88**, 1213 (1992).

[42] A. Tonnelier, H. Belmabrouk, and D. Martinez, *Neural Computation* **19**, 3226 (2007).

[43] D. Golomb, *Scholarpedia* **2**, 1347 (2007).

[44] E. Ullner and A. Politi, *Physical Review X* **6**, 011015 (2016).

[45] The coefficient of variation $cv(i)$ for the neuron $i$ is the ratio between the standard deviation and the mean of the interspike intervals associated with its firing activity. $CV$ is the ensemble average of the single neurons $cv(i)$.

[46] The value of $\overline{CV}$ is $\approx 0.8 - 0.9$ in the whole range of in-degrees $K \in [10:500]$ for $i_0/g_0^2 = 0.00055$ [45].

[47] Cluster activation is not related to clustering instabilities

observable in noisy globally coupled neural networks [52, 53].

[48] V. V. Klinshov and S. Y. Kirillov, *Physical Review E* **106**, L062302 (2022).

[49] V. Klinshov, P. Smelov, and S. Y. Kirillov, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **33** (2023).

[50] A. Torcini and P. Politi, *The European Physical Journal B-Condensed Matter and Complex Systems* **25**, 519 (2002).

[51] R. Albert and A.-L. Barabási, *Reviews of modern physics* **74**, 47 (2002).

[52] N. Brunel and D. Hansel, *Neural Computation* **18**, 1066 (2006).

[53] Y. Feld, A. K. Hartmann, and A. Torcini, *New Journal of Physics* **26**, 063017 (2024).