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

HAL
open science

# How Random Connectivity Shapes the Fluctuating Dynamics of Finite-Size Neural Populations

Nils E Greven, Jonas Ranft, Tilo Schwalger

► To cite this version:

Nils E Greven, Jonas Ranft, Tilo Schwalger. How Random Connectivity Shapes the Fluctuating Dynamics of Finite-Size Neural Populations. PRX Life, 2026, 4 (1), pp.013007. 10.1103/shvm-x4x6. hal-05513579

HAL Id: hal-05513579

https://cnrs.hal.science/hal-05513579v1

Submitted on 16 Feb 2026

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 CC BY 4.0 - Attribution - International License
# How Random Connectivity Shapes the Fluctuating Dynamics of Finite-Size Neural Populations

Nils E. Greven<sup>1,2</sup> Jonas Ranft<sup>3,*</sup> and Tilo Schwalger<sup>1,2,†</sup>

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

(Received 20 December 2024; accepted 5 December 2025; published 14 January 2026)

Mesoscopic models of finite-size neuronal populations are crucial to understand the dynamics of neural networks in the brain, especially their fluctuations and response to stimuli. However, current theories to derive such models are based on the assumption of homogeneous all-to-all (full) connectivity, neglecting the variance in the connectivity of biologically realistic networks with connection probabilities $p < 1$ (nonfull connectivity). To gain insight into the different fluctuation mechanisms underlying the neural variability of populations of spiking neurons, we derive and analyze a stochastic mean-field model for finite-size networks of Poisson neurons with random connectivity (including nonfull connectivity), external noise, and disordered mean inputs. We treat the quenched disorder of the connectivity by an annealed approximation enabling a doubly stochastic description of synaptic inputs for finite network size. A further reduction leads to a low-dimensional closed system of coupled Langevin equations for the mean and variance of the neuronal membrane potentials as well as a variable capturing finite-size fluctuations arising specifically in the case of connectivity disorder. Comparing to microscopic simulations, we find that the mesoscopic model describes the fluctuations and nonlinearities well and outperforms previous mesoscopic models that neglected the variance in the connectivity. The joint effect of connectivity disorder and finite network size can be analytically understood by a softening of the effective nonlinearity and the multiplicative character of spiking noise. The mesoscopic theory shows that quenched disorder can stabilize the asynchronous state, and it correctly predicts large quantitative and nontrivial qualitative effects of connection probability on the variance of the population firing rate and its dependence on stimulus strength. Our theory thus elucidates how disordered connectivity shapes nonlinear dynamics and fluctuations of neural populations at the mesoscopic scale and showcases a useful mean-field method to treat random connectivity in finite-size, spiking neural networks.

DOI: 10.1103/shvm-x4x6

## I. INTRODUCTION

Biological systems often exhibit significant fluctuations and variability in their dynamics in contexts ranging from gene expression to neural activity in the brain. The possible roles of fluctuations for the dynamics and functions of such systems have thus received increasing attention [1–7]. An important tool to understand the collective dynamics of biological systems consisting of many interacting units are low-dimensional descriptions of the effective dynamics at the mesoscopic or macroscopic scale. Taking fluctuations properly into account can be necessary to establish a correct understanding of the system and can provide crucial insight into the origins of observed biological variability. At the mesoscopic scale, fluctuations naturally arise from the finite number of underlying units. These “finite-size” fluctuations may then be accounted for by stochastic low-dimensional

models. In the context of living systems, examples range from Langevin equations for gene regulation [8] to coupled molecular motors [9–11] to ecology [12].

In neuroscience, neural-population models (also known as firing-rate or neural-mass models), such as the Wilson-Cowan model [13], have been highly successful in describing various cortical phenomena such as oscillations [14], multistability [15], and nonlinear response properties [16]. Stochastic versions of neural-population models have been critical for understanding fluctuating neural population dynamics such as metastability [15,17,18], stochastic oscillations [19], and response properties of cortical variability [20,21]. Despite their mathematical tractability, however, these models are often phenomenological, lacking a clear link to an underlying finite-size network of spiking neurons.

A link to biophysical properties is maintained through bottom-up mean-field approaches [22–28], where the macroscopic or mesoscopic dynamics is analytically derived from a microscopic model. The resulting mean-field dynamics are amenable to mathematical analysis [29–33] and enable efficient and accurate simulations of large spiking neural networks [34–37]. Recently, a bottom-up theory for the mesoscopic dynamics of finite-size populations of generalized integrate-and-fire neurons with escape noise has been developed [33,36]. It suggests an efficient multiscale modeling framework of cortical circuits [36,38,39] grounded in realistic

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote n="1"><sup>1</sup>Institute of Mathematics, Technical University Berlin, 10623 Berlin, Germany</footnote>

<footnote n="2"><sup>2</sup>Bernstein Center for Computational Neuroscience Berlin, 10115 Berlin, Germany</footnote>

<footnote n="3"><sup>3</sup>Institut de Biologie de l'ENS, Ecole normale supérieure, PSL University, CNRS, Paris, France</footnote>

<footnote>*Contact author: jonas.ranft@ens.psl.eu</footnote>

<footnote>†Contact author: tilo.schwalger@bccn-berlin.de</footnote>
neuron dynamics [40,41]. However, this and other mean-field theories [33, 36, 42–46] for finite-size neural networks typically assume homogeneous, all-to-all (full) connectivity. This assumption can be regarded as an approximation of a nonfully connected random network by its homogeneous mean connectivity while neglecting its quenched variability [36,47]. On the other hand, the role of random connectivity for intrinsically generated fluctuations of synaptic inputs has been extensively studied using dynamical mean-field theories [29, 30, 48–53], albeit only in the limit of infinitely large networks. Both assumptions—homogeneous full connectivity ignoring quenched variability and infinite network size ignoring fluctuations—are severe limitations for describing biologically realistic, recurrent neural networks. For example, rat barrel cortex contains only on the order of $N \sim 100$ to $N \sim 1000$ neurons per cell type per layer [54]—a population size for which finite-size fluctuations can be significant, as shown in computational models [36,55–59]. Furthermore, the local connection probabilities $p$ in mouse visual cortex vary on a wide range ($0 \lesssim p \lesssim 0.7$) with low average ($p \sim 0.1$) and larger connection probabilities ($p \sim 0.5$) for similarly tuned neurons [60]. How such pronounced deviations from fully connected networks ($p = 1$) affect the mesoscopic dynamics of finite-size neural populations is still poorly understood.

Both finite population size ($N < \infty$) and quenched random connectivity (e.g., $p < 1$) cause intrinsically generated (endogenous) fluctuations of the recurrent, synaptic input currents. These intrinsic fluctuations act as multiplicative noise whose strength increases with population firing rates. On the other hand, neural variability in the cortex often decreases upon stimulation [61]—a phenomenon that could be explained, among others, with a model where variability is generated by an external, additive noise [21]. To clarify the role of intrinsic fluctuations for such cortical phenomena under biologically realistic neuron numbers and connection probabilities requires a mean-field theory that accounts for these different origins of fluctuations.

In a homogeneous, fully connected network, as arising under a mean-connectivity approximation, the input fluctuations are common to all neurons. This is because all neurons receive identical recurrent inputs and thus experience the same finite-size fluctuations of the population activity. These “coherent” fluctuations lead to stochastic population equations with finite-size noise of order $1/\sqrt{N}$. In contrast, in random networks with connection probability $p < 1$, neurons no longer receive identical synaptic inputs, and hence input fluctuations become partially decorrelated. Consequently, the assumption of perfectly coherent fluctuations no longer holds, raising questions about the validity of the mean-connectivity approximation for finite-size networks.

A mean-field description of random networks (including nonfully connected networks, $p < 1$) requires the treatment of fluctuations resulting from a “quenched” (i.e., temporally constant) random synaptic connectivity. Previous approaches to random networks considered the macroscopic limit $N \to \infty$. Depending on the concrete limiting procedure and network architecture, the recurrent fluctuations seen by different neurons either vanish [22, 28, 31, 32] or become perfectly independent (“incoherent”) in this limit [29, 30, 48, 50–53] unless correlations are induced externally. A famous example is the

mean-field theory for sparsely connected networks by Brunel and coworkers [29,30]. In this theory, the neural population dynamics is deterministic and given by a nonlinear Fokker-Planck equation. The effect of incoherent fluctuations of the recurrent synaptic input caused by the random connectivity appears as an additional contribution to the diffusion coefficient in this equation. Although this approach relies on an annealing approximation that neglects temporal correlations in the spike trains caused by the quenched disorder, the nonlinear Fokker-Planck equation has been successful to capture the influence of intrinsic fluctuations on the nonlinear population dynamics (for a treatment of temporal correlations, see, e.g., Ref. [62]). Using this theory, intrinsic fluctuations have been shown to crucially shape nonlinear response properties [63], network oscillations [29,64], and multistability [56].

Even though dynamical mean-field theories for $N \to \infty$ have been used for real neural networks of finite size (see, e.g., Ref. [65]), they cannot describe finite-size fluctuations at the population level. Furthermore, finite-size networks are necessarily outside the sparse limit (i.e., $p > 0$) which entails a nonvanishing probability of shared recurrent inputs among neurons. Therefore the assumption of perfectly incoherent noise in dynamical mean-field theories cannot hold fully true for $N < \infty$. While finite-size extensions of the mean-field theory for sparsely connected networks have been proposed in Refs. [29,30,66], low-dimensional neural population models that reveal and take into account the distinct effects of coherent and incoherent fluctuations are currently lacking.

In this paper, we consider finite-size networks with random connectivity, allowing arbitrary connection probabilities $p \in [0, 1]$. For these networks, we derive and analyze a simple mesoscopic dynamics that describes incoherent and coherent recurrent fluctuations. To this end, we study a network of Poisson neurons using an annealed approximation [29,67] of the quenched random connectivity. The annealed approximation neglects temporal correlations of incoherent fluctuations but captures the effects of random connectivity and finite network size surprisingly well. In contrast to the case $N \to \infty$ of previous theories [29,30], we need to treat the fluctuating input as a doubly stochastic process to account for both finite-N and random connectivity. This description eventually yields a three-dimensional Langevin equation for the mesoscopic dynamics. We compare this dynamics with a naive mean-field theory in which the variance of the random connectivity is neglected, corresponding to the ad hoc mean-field approximation used in previous studies [36,47]. While the focus of the paper is the Langevin dynamics for a homogeneous population with quenched or annealed random connectivity and common external noise, the theoretical framework can be extended to heterogeneous populations, where the resting potentials or external currents of the neurons exhibit quenched Gaussian disorder. In summary, we present a simple theoretical framework for the effects of various sources of noise and quenched disorder on the neural population dynamics at the mesoscopic scale ($1 \ll N < \infty$).

The paper is organized as follows. We introduce the mesoscopic network model with quenched random connectivity and its annealed approximation in Sec. II. By reducing the annealed network to a mean-connectivity network with incoherent dynamical noise, we derive in Sec. III a mesoscopic
population model as a system of three coupled stochastic differential equations. In Sec. IV, we then analyze the mesoscopic model with respect to fixed-point solutions and their stability, the linear response to dynamic stimuli and second-order statistics. In particular, we compare the obtained analytical results to the corresponding statistics of the quenched and annealed microscopic models using extensive simulations. We conclude and discuss our results in Sec. V. The extension to quenched Gaussian disorder as well as longer, detailed calculations are provided in the Appendix.

## II. MICROSCOPIC MODELS

### A. Quenched network model

At the microscopic level, we consider a random network of N interacting Poisson neurons (also called nonlinear Hawkes processes). While our theory will be developed for general random connectivity with finite variance, we focus on the specific example of binary, nonfull connectivity to study the effect of connection probabilities $p < 1$. The Poisson neurons fire spikes stochastically with conditional intensity $\lambda_i(t) = \phi(h_i(t^{-})), i = 1, \dots, N$, where $h_i(t)$ is the input potential of neuron i and $\phi$ is a non-negative hazard function. The dynamics of the input potentials is given by a system of coupled first-order equations with delay,

$$\tau \frac{dh_i}{dt} = -h_i + \mu(t) + \frac{w}{C} \sum_{j=1}^{N} a_{ij} \dot{Z}_j(t-d), \quad (1)$$

$i = 1, \dots, N$. Here $\tau$ is the time constant of the low-pass filter dynamics, $\mu(t)$ is an external drive, and $\dot{Z}_i(t) = \sum_{j=1}^{N} \delta(t - t_{i,l})$ is the spike train of neuron i with $\{t_{i,l}\}$ being the individual spike times. The presence or absence of synaptic connections between neurons are described by the random adjacency matrix $[a_{ij}]$ with fixed in-degree $C = \sum_{j=1}^{N} a_{ij}$ and elements $a_{ij} \in \{0, 1\}$ that are marginally Bernoulli distributed with mean $\langle a_{ij} \rangle = p = C/N$ representing the connection probability. This random connectivity can be constructed by choosing for each neuron C presynaptic neurons randomly, and the resulting connectivity is fixed (“quenched”) in time. Furthermore, $d$ is the synaptic delay and $w = CJ$ is the total coupling strength with $J$ being the efficacy of a single synapse (i.e., $J/\tau$ is the jump size in millivolt of the postsynaptic potential caused by a single presynaptic spike). Since we focus on inhibitory networks in this paper, we assume that $J < 0$.

As a concrete hazard function, we choose here a sigmoidal function in the form of an error function [22]:

$$\phi(h) = r_m \Phi(\beta(h - \vartheta)), \quad (2)$$

where $\Phi(x) = [1 + \erf(x/\sqrt{2})]/2$. Although not essential for the general theory, this choice will allow us to analytically calculate the first and second moments of $\phi(h_i)$. The parameter $\beta$ determines the steepness around the inflection point. Without loss of generality, we choose the position of the inflection point $\vartheta$ to be zero because we can always measure voltages with respect to the potential $\vartheta$. The sigmoidal shape with an upper limit $r_m$ prevents the Poisson neurons to fire with arbitrarily high rates. While the concave behavior and saturation of the sigmoidal hazard function for $h > 0$ is technically useful, we are in the following mainly interested in dynamical

regimes operating in the convex part of the hazard function $\phi(h)$, i.e., corresponding to mean inputs below the inflection point, $h < 0$. This is because biologically realistic hazard functions of cortical neurons are typically convex [68,69]. As default parameters of the model, we choose $\tau = 20$ ms, $r_m = 100$ Hz and $d = 0$ ms unless otherwise noted.

The aim of this paper is to derive a mesoscopic population model that generates population activities that statistically match the population activities obtained from a microscopic network simulation. We define the population activity $A_N(t, \Delta t)$ with respect to a time discretization with time step $\Delta t$ as the total number of spikes per neuron and time step:

$$A_N(t; \Delta t) = \frac{\Delta Z(t)}{N \Delta t}. \quad (3)$$

Here $\Delta Z(t)$ denotes the total number of spikes in the time interval $(t, t + \Delta t]$.

### B. Mean-connectivity network

A major obstacle in deriving population models analytically is the quenched randomness of the connectivity $J_{ij} := J_{a_{ij}}$. A standard mean-field approach to tackle this problem is to neglect fluctuations in the recurrent synaptic input by scaling the synaptic weights inversely proportional to the network size, $J_{ij} = w_{ij}/N$, where the coupling strength $w_{ij}$ is assumed to be of order 1. Then, in the limit $N \to \infty$, the synaptic input $N^{-1} \sum_j w_{ij} \dot{Z}_j$ converges to the deterministic input $\langle w_{ij} \rangle \langle \dot{Z}_j \rangle$ [27,28], where $\langle \cdot \rangle$ denotes the average over both quenched disorder and Poisson noise. In our case, this corresponds to a network where the quenched connectivity $w_{ij}$ [Fig. 1(a), left] is replaced by the mean connectivity $\langle w_{ij} \rangle = NJ \langle a_{ij} \rangle = CJ = w$ [Fig. 1(b), left]. Importantly, the mean connectivity effectively yields a fully connected network with homogeneous coupling $w$. In the limit $N \to \infty$, the synaptic input is thus not affected by the variance

$$\sigma_w^2 = w^2(1-p)/p \quad (4)$$

of the disordered couplings $w_{ij}$ but only by the mean $\langle w_{ij} \rangle$.

The replacement of all elements of the adjacency matrix $a_{ij}$ with their mean $p$ has been used ad hoc also for finite network size $N < \infty$ [36,47]. However, this ad hoc approach is no longer exact for finite $N$ but must be regarded as an approximation that retains the mean synaptic connectivity $\langle J_{ij} \rangle = pJ$ but ignores all higher-order cumulants of the random variables $J_{ij}$. In the following, we will refer to this effective, fully connected network as the mean-connectivity network, whereas the original model, Eq. (1), will be referred to as the quenched network (Table I). Replacing the connection strengths $w_{ij}$ with their means $w$, we obtain from Eq. (1) the dynamics of the mean-connectivity network:

$$\tau dh_i(t) = [-h_i + \mu(t)]dt + \frac{w}{N} \sum_{j=1}^{N} dZ_j(t-d), \quad (5)$$

where $dZ_j(t) \sim \text{Pois}(\phi(h_j(t))dt)$ are the conditionally independent Poisson increments of the spike count of neuron j. The mean-connectivity network of Poisson neurons, Eq. (5), is a microscopic model that admits an exact reduction to a mesoscopic model for the total spike count $dZ(t) = \sum_{i=1}^{N} dZ_i(t)$.
<figure><img src="image_4.png" /><figcaption>FIG. 1. Response of the population rate variance to a change of the mean stimulus strength in the original, quenched random network (a) and in a corresponding fully connected, mean-connectivity network that only accounts for the mean synaptic weight (b). Top: Noisy external stimulus $\mu(t)$ whose mean increases at $t = 0.08$ s in a steplike fashion from $\mu_0 = 28$ mV to $\mu_0 = 50$ mV; the noise strength $\sigma_{\text{ext}}^2 = 1$ mV$^2$ is kept fixed. (a) Schematic of a quenched network with fixed adjacency matrix $a_{ij}$ (left) and the trial-averaged variance of the population firing rate as a function of time (right) for different network sizes (green line, $N = 5 \times 10^4$; blue line, $N = 10^3$; in both cases $C = 100$). (b) Same for a mean-connectivity network, Eq. (5), as an ad hoc approximation of the quenched network, where the quenched connectivity $a_{ij}$ is replaced by the average connectivity $\langle a_{ij} \rangle$. While the mean stationary population rates after the step are equal in both models (measured rates are provided in the right panels), the variances of the population rate strongly differ (note the different scaling of the y axis). Other parameters: $\beta = 5$ mV$^{-1}$, $w = -1$ mV s.</figcaption></figure>

To see this, we note that the sum of independent Poisson numbers in Eq. (5) is a Poisson number with mean $Nr(t-d)dt$, where

$$r(t) = \frac{1}{N} \sum_{i=1}^{N} \phi(h_i(t)) \quad (6)$$

<table><caption>TABLE I. Summary of microscopic network models.</caption><thead><tr><th>Network</th><th>Equation</th><th>Connectivity</th></tr></thead><tbody><tr><td>Quenched</td><td>(1)</td><td>Random, fixed in time,<br>⟨w<sub>ij</sub>⟩ = w, var(w<sub>ij</sub>) = σ<sub>w</sub><sup>2</sup></td></tr><tr><td>Mean-connectivity</td><td>(5)</td><td>Mean-matched all-to-all,<br>⟨w<sub>ij</sub>⟩ = w, var(w<sub>ij</sub>) = 0</td></tr><tr><td>Annealed</td><td>(8), (9)</td><td>Random, dynamic<br>(resampling in time),<br>⟨w<sub>ij</sub>(t)⟩ = w, var(w<sub>ij</sub>(t)) = σ<sub>w</sub><sup>2</sup></td></tr></tbody></table>

is the stochastic population rate. Furthermore, because all neurons are driven by the same synaptic input, their input potentials $h_i(t)$ coincide (or converge after an initial transient). The common input potential $\tilde{h}(t)$ thus obeys the mesoscopic dynamics

$$\tau d\tilde{h}(t) = [-\tilde{h} + \mu(t)]dt + \frac{w}{N} dZ(t - d), \quad (7)$$

where now $dZ(t) \sim \text{Pois}(N\phi(\tilde{h}(t))dt)$ describes the increment of the total spike count of the population. Equation (7) describes the mesoscopic dynamics as it only involves the population activity $A_N(t, \Delta t) = \int_t^{t+\Delta t} dZ(t)/(N\Delta t)$.

Simulations of the quenched and the mean-connectivity network yield similar results for the stationary mean firing rates and mean input potentials. However, the mean-connectivity network may differ drastically compared to the quenched network when it comes to the nonstationary response (linear response function) and the fluctuations (second-order statistics). To illustrate the failure of the mean-connectivity approximation, let us consider the variance of the stochastic population rate $r(t)$. This variance represents the part of the mesoscopic neural variability [variance of the population activity $A_N(t, \Delta t)$] that is caused by the rate variability rather than the Poisson spiking noise. It can be measured from population activity data (see Appendix C) or directly computed from a simulated time series of $r(t)$ in our model. A nonvanishing variance of the population rate can occur for two reasons in our model: first, because of intrinsic finite-size fluctuations and, second, because of externally injected common noise. Here we model the common external noise as a Gaussian white noise [more precisely, $\mu(t) = \bar{\mu}(t) + \sqrt{\tau\sigma_{\text{ext}}^2}\hat{\zeta}(t)$, where $\bar{\mu}(t)$ is the mean stimulus and $\hat{\zeta}(t)$ is a standard Gaussian white noise process].

Interestingly, simulations of the mean-connectivity network [Eqs. (5) or (7)] strongly overestimate the population-rate variance of the quenched network by more than one order of magnitude (Fig. 1, note the different scale of the y axes). We stress that this difference is not an effect of a difference in firing rates: The stationary mean population rates $\langle r \rangle$ are the same in both networks [e.g., the numerical values of the firing rates at high stimulus amplitudes are approximately equal as indicated in Figs. 1(a) and 1(b)]. In addition to this large quantitative difference, we also observe a qualitatively different response to fast stimulus changes when the network is large ($N = 5 \times 10^4$, green lines in Fig. 1): When the mean stimulus strength jumps from 28 to 50 mV, the variance of the population rate exhibits a small but significant decrease in the quenched network [Fig. 1(a), right]. This suppression of variability is in marked contrast to the prediction by the mean-connectivity network, which exhibits instead an increase in the variance [Fig. 1(b), right]. This qualitatively different behavior of the population-rate variance in response to a step stimulus is, however, not observed for smaller network size ($N = 10^3$, but otherwise same parameters and similar firing rates; blue lines in Fig. 1).

How can we understand nonstationary responses and neutral variability of mesoscopic variables theoretically? In this paper, we present a “second-order” mesoscopic mean-field theory that also captures second-order statistics and non-stationary responses quantitatively and explains nontrivial,
<figure><img src="image_5.png" /><figcaption>FIG. 2. Schematic illustration of the annealed network and its approximation by a mean-connectivity network with additional noise. (a) In the annealed network model, Eq. (8), the random connectivity is resampled at each spike time resulting in an adjacency matrix $a_{ij}$ that varies at each spike $k$. (b) The annealed model can be mapped to a mean-connectivity network with additional noise, Eq. (12). The network is all-to-all connected with an average connectivity $\langle a_{ij} \rangle$. Additionally, each neuron is subject to an independent noise that captures the random connectivity. (c) Example trajectories of the input potentials in the mean-connectivity network with additional noise (gray). The distribution $P$ of the input potentials follows a Gaussian distribution with a time-dependent mean $\bar{h}(t)$ and a time-dependent variance $\sigma^2(t)$. The second-order mean-field model (20) aims to describe the stochastic dynamics of these two mean-field variables.</figcaption></figure>

qualitative phenomena such as the suppression of variability observed in Fig. 1(a). In particular, the explanation of the failure of the mean-connectivity approximation with respect to the stationary variances in Fig. 1(b) and its resolution by the “second-order” mean-field theory will be presented towards the end of the paper in Sec. IV E.

### C. Annealed network model

The quenched connectivity causes heterogeneity of neural firing with spatial and temporal correlations which often make a direct mean-field treatment infeasible. Here we follow a different approach based on a dynamical resampling of the connectivity (“annealing”) at each spike [29,67] [Fig. 2(a)]. This procedure largely ignores the part of the spatial and temporal correlations that is caused by the quenched connectivity but leads to a tractable model of a fully connected neural network with additional dynamical noise that, as we will see, retains and captures the main statistical features of the recurrent synaptic input to a neuron [Fig. 2(b)]. We expect that the annealed approximation is better in dense networks of strongly correlated neurons. In this case, the temporal correlations in the spike input of a given neuron are partially retained upon randomization of the presynaptic neurons. In contrast, in sparse networks, the membrane potentials are independent among neurons, and therefore, the resampling of presynaptic neurons in time destroys the temporal correlations in the presynaptic spike train of a fixed presynaptic neuron.

In the annealed approximation, the connections $a_{ij}$ are resampled at each spike occurring in the network. Thus, the time-varying adjacency matrix can be written as $[\hat{a}_{ijk}]$, where the additional index $k$ labels the $k$th spike occurring in the network. Then, for a fixed $k$, the $N \times N$-adjacency matrix $[\hat{a}_{ij}]_k$ is still a random matrix with fixed in-degree and therefore has independent rows but dependent elements $\hat{a}_{ijk}$ within a given row $i$. However, the dependence caused by the fixed in-degree constraint is immaterial for the recurrent inputs in the annealed approximation. This is because only the $j_k$th column of the matrix $[\hat{a}_{ij}]_k$, which has independent elements, matters for the transmission of the $k$th spike, where $j_k$ is the index of the neuron that fired the $k$th spike, and because of the resampling of the adjacency matrix. The latter implies that, for different spike times, the columns are independent. From the perspective of the recurrent inputs, we can thus equally well resample the elements of the adjacency matrix as independent Bernoulli variables $a_{ijk} \sim \text{Ber}(p)$ with success probability $p = C/N$. Thus, our model, Eq. (1), changes in the annealed approximation to

$$ \tau dh_i(t) = [-h_i + \mu(t)]dt + \frac{1}{N}dY_i(t-d), \quad (8) $$

$$ \frac{dY_i}{dt}(t) = \frac{w}{p} \sum_{j=1}^{N} \sum_{l} a_{ijk_{j,l}} \delta(t - t_{j,l}), \quad (9) $$

where the index $k_{j,l}$ delivers the location of the $l$th spike time $t_{j,l}$ of neuron $j$ within the sequence of all spikes in the network. The quantity $Y_i(t)$ is the integrated synaptic input received by neuron $i$ until time $t$. We refer to this model as the annealed network (Table I). We note that in simulations of the annealed network, it is not necessary to resample the entire adjacency matrix at each spike but only the column that corresponds to the neuron that fired the spike. In other words, whenever a neuron fires a spike, this spike is transmitted independently to postsynaptic neurons with probability $p = C/N$.

By generalizing from binary weights $w_{ij} = \frac{w}{p}a_{ij}$ to an arbitrary connectivity matrix $w_{ij}$, the annealed network can also be defined for such general connectivity through dynamic resampling, which produces time-dependent synaptic weights $w_{ij}(t)$. In the specific case of binary synapses, $w_{ij}(t_{j,l}) = \frac{w}{p}a_{ijk_{j,l}}$ as in Eq. (9).

Apart from the random connectivity, another source of quenched disorder that can be included in our theory is heterogeneity in the external current $\mu_0$. For the sake of clarity, we will begin our analysis without heterogeneity in $\mu_0$ and will defer the straightforward extension to heterogeneity to a brief section in Appendix A. The parameter values of the microscopic model used in simulations are provided in Table II

## III. DERIVATION OF THE MESOSCOPIC MODEL

### A. A mean-connectivity network with additional noise that accounts for second-order connectivity statistics

The first step of the derivation is a temporal coarse-graining of the annealed model leading to an effective, fully connected network with homogeneous weights (mean-connectivity network) and additional noise. The additional noise will be the
crucial difference to the mean-connectivity network discussed in Fig. 1. To coarse-grain, we integrate Eq. (8) over a small time step of length $\Delta t$:

$$\tau \Delta h_i(t) = \int_t^{t+\Delta t} [-h_i(s) + \mu(s)] ds + \frac{1}{N} \Delta Y_i(t-d). \quad (10)$$

Here $\Delta h_i(t) = h_i(t + \Delta t) - h_i(t)$ is the increment of the input potential, $\Delta Y_i(t) = \sum_{j=1}^{N} w_{ij}(t) \Delta Z_j(t - d)$ is the synaptic input received by neuron $i$ during the time step $\Delta t$ and $\Delta Z_j(t) = 1$ if neuron $j$ has a spike during that time step and $\Delta Z_j(t) = 0$ otherwise. In the annealed approximation, the time-varying coupling strengths $w_{ij}(t)$ are independent random variables (across time steps and synapses) with mean $w$ and variance $\sigma_w^2$. In our specific example, the random variables are binary taking the value $w/p$ with probability $p$ and $0$ otherwise. However, the derivation does not hinge upon binary connectivity but applies more generally to any distribution of coupling strengths with finite variance $\sigma_w^2 < \infty$.

The synaptic input $\Delta Y_i(t)$ can be regarded as the result of a doubly stochastic process: On the one hand, given the total spike count $\Delta Z(t) = \sum_{i=1}^{N} \Delta Z_i(t)$, each $\Delta Y_i(t)$ is the sum of $\Delta Z(t)$ independent random variables. Hence, the conditional mean and covariance are

$$\langle \Delta Y_i(t) | \Delta Z(t) \rangle = w \Delta Z(t), \\ \text{cov}(\Delta Y_i(t), \Delta Y_j(t) | \Delta Z(t)) = \sigma_w^2 \Delta Z(t) \delta_{ij},$$

where $\delta_{ij}$ denotes the Kronecker delta. If $N$ is large such that $p\Delta Z(t) \gg 1$, then the synaptic input $\Delta Y_i(t)$ for a given $\Delta Z(t)$ is approximately Gaussian distributed. Under this condition, we can write

$$\Delta Y_i(t) \approx w \Delta Z(t) + \sigma_w \sqrt{\Delta Z(t)} n_i(t), \quad (11)$$

where $n_i(t) \sim N(0, 1)$ are independent standard normal random variables. The first term in Eq. (11) is proportional to the total number of spikes or the mesoscopic population activity and thus corresponds to the input in an effective mean-connectivity network. This term is common to all neurons. In contrast, the second term is independent across neurons and only arises when $\sigma_w > 0$. Specifically, it captures the individual differences of the synaptic inputs in a nonfully connected (“diluted”) network with $p < 1$.

On the other hand, $\Delta Z(t)$ is itself a Poisson random variable with mean $Nr(t)\Delta t$. From Eq. (11), or, alternatively, from the law of total expectation and total covariance, we find that the synaptic inputs $\Delta Y_i(t)$ exhibit the following mean and covariance for a given stochastic population rate $r(t)$:

$$\langle \Delta Y_i(t) | r(t) \rangle = w Nr(t) \Delta t \\ \text{cov}(\Delta Y_i(t), \Delta Y_j(t) | r(t)) = (w^2 + \sigma_w^2 \delta_{ij}) Nr(t) \Delta t.$$

The last equation shows that pairwise correlations of the synaptic input to different neurons are nonzero for finite $N$. These correlations arise from the common fluctuations of the first term in Eq. (11) representing the population activity.

The covariance structure of the synaptic inputs $\Delta Y_i(t)$ is preserved if their representation by Eq. (11) is replaced by

$$\Delta \hat{Y}_i(t) = w \Delta Z(t) + \sigma_w \sqrt{Nr(t) \Delta t} n_i(t).$$

This representation allows to take the temporal continuum limit. Substituting $\Delta \hat{Y}_i(t-d)$ for $\Delta Y_i(t-d)$ in Eq. (10), dividing by $\Delta t$ and taking the limit $\Delta t \to 0$, we obtain

$$\tau \frac{d h_i}{d t} = -h_i + \mu(t) + \frac{w}{N} \sum_{j=1}^{N} \dot{Z}_j(t-d) + \sigma_w \sqrt{\frac{r(t-d)}{N}} \zeta_i(t). \quad (12)$$

Here we have used that the population activity, Eq. (3), can be rewritten as

$$A_N(t, \Delta t) = \frac{1}{N \Delta t} \int_t^{t+\Delta t} \sum_{i=1}^{N} dZ_i(t'), \quad (13)$$

where the spike trains $\dot{Z}_i(t)$ have stochastic intensities $\phi(h_i(t^{-}))$. Furthermore, $\zeta_i(t)$ are independent Gaussian white noises obeying $\langle \zeta_i(t) \zeta_j(t') \rangle = \delta_{ij} \delta(t - t')$.

Equation (12) represents an effective, fully connected network with rescaled synaptic efficacy $w/N = pJ$ and additional dynamic noise whose intensity is proportional to the variance $\sigma_w^2$ of the coupling strengths $w_{ij}$. Therefore, the dynamics is the same as the mean-connectivity network, Eq. (5), but with an additional noise term that captures the fluctuations caused by the random connectivity. This representation also reveals the distinct effects of finite network size ($N < \infty$) and random dilution ($p < 1$) on the recurrent synaptic fluctuations. First, the stochastic input spikes from all neurons in the network contribute a common synaptic input current proportional to the population activity [third term on the right-hand side of Eq. (12)], which is stochastic for finite $N$ and thus exposes the coherent part of the synaptic fluctuations. More precisely, the population activity has conditional mean $r(t)$ [given all input potentials $h_i(t^{-})$, $i = 1, \dots, N$] and finite-size fluctuations of order $N^{-1/2}$ because the spike trains $\dot{Z}_i(t)$ are conditionally Poisson [given $h_i(t^{-})$] with conditional mean $r(t)$.

In contrast, the last term of Eq. (12) contains a Gaussian white noise $\zeta_i(t)$ that is independent for each neuron, and thus represents the incoherent part of the synaptic fluctuations. This incoherent noise captures the effect of the random connectivity ($\sigma_w^2 > 0$) and, specifically, the random dilution ($p < 1$). It therefore vanishes in the limit of fully connected, homogeneous networks ($p \to 1$, $\sigma_w^2 \to 0$), in which case the mean-connectivity network is recovered.

In the following, we assume a large network size $N$ such that a diffusion approximation of the shot noise can be made: Noting that $\Delta Z(t)$ is conditionally Poisson with mean and variance equal to $Nr(t)\Delta t$, the population activity can be approximated for large $N$ as

$$A_N(t; \Delta t) \approx \frac{\int_t^{t+\Delta t} A(t') dt'}{\Delta t}, \quad A(t) = r(t) + \sqrt{\frac{r(t)}{N}} \eta(t), \quad (14)$$

where $\eta(t)$ is a common Gaussian white noise with $\langle \eta(t) \eta(t') \rangle = \delta(t - t')$. In this form, the finite-size fluctuations of order $1/\sqrt{N}$ become explicit. We note that the diffusion approximation of the population activity, Eq. (14), can in principle become negative. However, if $N$ is sufficiently
large such that $Nr(t)\Delta t \gg 1$, then the probability that such event happens is negligibly small.

## B. Second-order mesoscopic mean-field dynamics

All-to-all uniformly connected networks with independent noise, such as the mean-connectivity network with additional noise derived in the previous section, have been well studied. In particular, in the limit $N \to \infty$, the population dynamics can be described by a Fokker-Planck equation [70], while for finite $N$, the stationary statistics of the fluctuations has been calculated using linear response theory [71]. However, despite earlier theoretical approaches [29,66], a derivation of stochastic population dynamics that captures both finite-size fluctuations and nonstationary dynamics remains challenging. In our setting, this derivation is facilitated by the absence of resets in the Poisson neuron model. To obtain a self-consistent description, we follow a standard two-step procedure: First, we consider the dynamics of the input potentials

$$ \tau \frac{d h_i}{dt} = -h_i + f(t) + g(t) \zeta_i(t), \quad (15) $$

driven by fixed, time-dependent input signals $f(t)$ and $g(t)$, and derive the dynamics of the mean and variance conditioned on these inputs. We then derive the dynamics of the stochastic rate $r(t)$ for fixed input signals $f$ and $g$ such that it depends on these input signals only through the conditional mean and variance and an additional colored noise. In the second step, in turn, the input signals will be expressed in terms of the stochastic rate $r(t)$ according to Eqs. (12), (13), and (14):

$$ \begin{aligned} f(t) &= \mu(t) + w \left[ r(t-d) + \sqrt{\frac{r(t-d)}{N}} \eta(t-d) \right], \\ g(t) &= \sigma_w \sqrt{\frac{r(t-d)}{N}}. \end{aligned} \quad (16) $$

This step will close the system for the conditional mean and variance and the colored noise because the stochastic rate depends itself only on these variables. In this way, we will obtain a self-consistent mean-field description for the annealed model, Eq. (12), as detailed in the following. In view of Eq. (16), we also note that the conditioning on the input signals $f$ and $g$ will allow for the interpretation of the mean-field variables at time $t$ as conditional averages given the history of the noise $\eta(t' - d)$, the stochastic rate $r(t' - d)$ and the external input $\mu(t')$ for $t' < t$.

In the first step, when $f$ and $g$ are fixed, the input potentials $h_i$ represent $N$ independent, time-inhomogeneous Ornstein-Uhlenbeck processes for which the ensemble mean $\bar{h}$ and ensemble variance $\sigma^2$ obey the dynamics

$$ \tau \frac{d \bar{h}}{dt} = -\bar{h} + f(t), \quad \tau \frac{d \sigma^2}{dt} = -2\sigma^2 + \frac{g^2(t)}{\tau}. \quad (17) $$

These variables are the conditional mean and variance described above. To relate the stochastic rate to these variables, we first consider the limit $N \to \infty$. In this limit, the stochastic

rate, Eq. (6), converges to its conditional mean

$$ \begin{aligned} F(\bar{h}(t), \sigma^2(t)) &:= \langle \phi(h_i(t)) | \bar{h}(t), \sigma^2(t) \rangle \\ &= \int_{-\infty}^{\infty} \phi(h) g_{\bar{h}, \sigma^2}(h) dh, \end{aligned} \quad (18) $$

where $g_{\bar{h}, \sigma^2}(h) = \exp[-(h - \bar{h})^2/(2\sigma^2)]/\sqrt{2\pi\sigma^2}$. Note that the conditional mean depends on the functions $f$ and $g$ only through the conditional mean $\bar{h}(t)$ and conditional variance $\sigma^2(t)$. Next, we consider the case of finite $N$. In this case, the stochastic rate $r(t)$ has finite-size fluctuations of order $1/\sqrt{N}$, hence we rewrite Eq. (6) as

$$ r(t) = \left[ F(\bar{h}(t), \sigma^2(t)) + \frac{1}{\sqrt{N}} \xi(t) \right]_+ \quad (19) $$

with $\xi(t) = \sqrt{N}[\frac{1}{N} \sum_{i=1}^N \phi(h_i(t)) - F(\bar{h}(t), \sigma^2(t))]$. The rectification bracket $[\cdot]_+$ = $\max(0, \cdot)$ could be inserted because $N^{-1} \sum_i \phi(h_i)$ is always a non-negative quantity. Including the rectification is necessary for a Gaussian approximation of $\xi$ to enforce a non-negative population rate $r(t)$, albeit negative values are extremely rare events for biologically relevant network sizes $N$. For fixed input functions $f$ and $g$, the second term $F(\bar{h}(t), \sigma^2(t))$ is deterministic, and hence, for large $N$, the variable $\xi(t)$ is an approximately Gaussian, colored noise with mean 0 and autocovariance function $c_{\bar{\xi}}^{f,g}(t, t') \approx \frac{c_{\bar{\xi}}^{f,g}(t,t')}{\sigma^2} c_{\bar{h}}^{f,g}(t, t')$, where $c_{\bar{h}}^{f,g}$ denotes the autocorrelation function of the process Eq. (15). Because $\xi(t)$ is Gaussian and $h_i(t)$ has an exponential autocorrelation function with correlation time $\tau$, we hence model the colored noise as an Ornstein-Uhlenbeck process of the form

$$ \tau \frac{d\bar{\xi}}{dt} = -\bar{\xi} + \sqrt{2\tau\sigma_{\lambda}^2(t)}\xi(t). $$

Here $\xi(t)$ is Gaussian white noise with $\langle \xi(t) \rangle = 0$ and $\langle \xi(t) \xi(t') \rangle = \delta(t - t')$, and

$$ \begin{aligned} \sigma_{\lambda}^{2}(t) &= \langle \phi(h_{i}(t))^{2} | \bar{h}(t), \sigma^{2}(t) \rangle - \langle \phi(h_{i}(t)) | \bar{h}(t), \sigma^{2}(t) \rangle^{2} \\ &= G(\bar{h}(t), \sigma^{2}(t)) \end{aligned} $$

is the conditional variance of the stochastic intensity $\lambda(t) = \phi(h_i(t))$ given by

$$ G(\bar{h}, \sigma^2) := \int_{-\infty}^{\infty} \phi^2(h) g_{\bar{h}, \sigma^2}(h) dh - F^2(\bar{h}, \sigma^2). $$

We note again that the conditional variance depends on the functions $f$ and $g$ only through the conditional mean $\bar{h}(t)$ and conditional variance $\sigma^2(t)$.

We now proceed with the second step of the mean-field derivation, where we close the system in the mean-field variables $\bar{h}$, $\sigma^2$, and $\xi$ self-consistently. To this end, the functions
f and g are related back to the stochastic mean-field rate $r(t)$ using Eq. (16). Inserting this equation into Eq. (17) we finally arrive at the mesoscopic mean-field (MF) dynamics

$$ 
\begin{aligned}
\tau \frac{d\bar{h}}{dt} &= -\bar{h} + \mu(t) + w \left[ r(t-d) + \sqrt{\frac{r(t-d)}{N}} \eta(t) \right] \\
\tau \frac{d\sigma^2}{dt} &= -2\sigma^2 + \frac{\sigma_w^2}{\tau N} r(t-d) \\
\tau \frac{d\xi}{dt} &= -\xi + \sqrt{2\tau G(\bar{h}(t), \sigma^2(t))\zeta(t)},
\end{aligned}
\quad (20) $$

where $r(t) = [F(\bar{h}(t), \sigma^2(t)) + \xi(t)/\sqrt{N}]_+$. The population activity $A_N(t, \Delta t)$ associated with this mesoscopic dynamics is given by Eq. (14). Because we consider the population mean $\bar{h}(t)$ and population variance $\sigma^2(t)$, we call Eq. (20) the second-order mean-field (MF2) dynamics.

The MF2 dynamics captures the effects of both finite $N$ and random connectivity. First, the finite-size effect manifests itself in the multiplicative noise term proportional to $\sqrt{r/N}$ in the equation for $\bar{h}$. The finite-size noise increases with larger firing rates, reflecting the Poisson noise in the underlying microscopic model, and it vanishes in the large-$N$ limit. Second, the effect of random connectivity becomes manifest in a nonvanishing variance $\sigma^2$, which is driven by the term proportional to $\sigma_w^2/N$. Note that in our case of binary random connectivity, the variance of the synaptic weights is $\sigma_w^2 = w^2(1-p)/p$ [cf. Eq. (4)]. Thus, for full connectivity, $p=1$, the variance $\sigma^2$ vanishes. It also vanishes for $N \to \infty$ unless $C$ is kept constant (sparse limit) in which case $\sigma_w^2/N = w^2(1/C - 1/N)$ remains of order 1. A nonvanishing variance $\sigma^2$ (in our case $p < 1$), in turn, impacts the population rate $r$ through the averaged nonlinearity $F$, Eq. (18). However, the population average, Eq. (6), is not perfectly equal to $F$ for finite $N$ and $\sigma^2 > 0$. The imperfect averaging leads to finite-size fluctuations $\xi(t)$—a third effect caused by the combination of both finite network size and random connectivity. This effect vanishes when $N \to \infty$ or $\sigma_w^2 = 0$ ($p=1$) because the noise strength is proportional to $\sqrt{G/N}$ and $G$ vanishes for $\sigma^2 = 0$.

If we had neglected the dispersion of the input potentials due to the random connectivity, i.e., assuming $\sigma^2 \equiv 0$, then we would have arrived at a first-order MF theory (MF1) with only one equation for the mean $\bar{h}(t)$:

$$ 
\begin{aligned}
\tau \frac{d\bar{h}}{dt} &= -\bar{h} + \mu(t) + w \left[ r(t-d) + \sqrt{\frac{r(t-d)}{N}} \eta(t) \right] \\
r(t) &= \phi(\bar{h}(t)).
\end{aligned}
\quad (21) $$

Here we have used that, for $\sigma = 0$, the population transfer function $F(\bar{h}, \sigma^2)$ simplifies to the transfer function of a single Poisson neuron $\phi(\bar{h}) = F(\bar{h}(t), 0)$ and that $G(\bar{h}, 0) = 0$. Theoretical predictions of MF1, Eq. (21), can always be obtained from the corresponding result of MF2 by taking the limit $p \to 1$ while keeping the network size $N = C/p$ and the coupling strength $w = JpN$ constant.

While the MF1 dynamics, Eq. (21), still captures finite-size effects through the noise term, it does not describe the effect of randomness in the connectivity $\sigma_w^2 > 0$ when $p < 1$. In fact, MF1 precisely corresponds to the mean-connectivity

approximation that underlies the assumption of full connectivity made in previous finite-size theories [33,36,42–47]. This can be seen by applying the diffusion approximation of the Poissonian shot noise for large $N$, Eq. (14), to the exact mesoscopic dynamics, Eq. (7), of the mean-connectivity network, Eq. (5), which recovers MF1. Thus, by comparing simulations of the microscopic models (quenched and annealed networks) with MF1 and MF2, we will thus be able to judge the mean-connectivity approximation used in previous studies [36,47] and the correction for the effect of random connectivity ($p < 1$) given by the second-order MF theory.

In the limit $N \to \infty$, the MF2 theory, Eq. (20), converges to either a one- or a two-dimensional, deterministic dynamics depending on whether the connectivity is dense or sparse. In the dense limit, where $p$ is kept constant and, hence, $C = pN \to \infty$, the variance $\sigma^2$ as well as the fluctuations $N^{-1/2}\eta$ and $N^{-1/2}\xi$ vanish, leaving a one-dimensional deterministic firing-rate model. This case is equivalent to the large-$N$ limit of the first-order MF theory and can thus neither capture finite-size nor $p < 1$ effects. In contrast, in the sparse limit, where $C$ is kept constant and, hence, $p = C/N \to 0$, the noise terms $N^{-1/2}\eta$ and $N^{-1/2}\xi$ vanish but the variance $\sigma^2$ does not. The result is a two-dimensional, deterministic mean-field dynamics for the mean and variance of the input potentials,

$$ 
\begin{aligned}
\tau \frac{d\bar{h}}{dt} &= -\bar{h} + \mu(t) + wr(t-d) \\
\tau \frac{d\sigma^2}{dt} &= -2\sigma^2 + \frac{w^2}{\tau C} r(t-d),
\end{aligned}
\quad (22) $$

where $r(t) = F(\bar{h}(t), \sigma^2(t))$. This sparse connectivity limit is similar to the classical model of Amari [22]. While it captures the effect of random connectivity, it cannot describe finite-size effects and hence second-order statistics caused by finite-size fluctuations. In the following, we will refer to the dynamics in the sparse limit, Eq. (22), in short as MFsparselim.

The effects described by the three population dynamics MF1, MF2, and MFsparselim for different parameter regimes are summarized in Table III (Appendix J).

### Functions F and G for an error-function nonlinearity

For the error function nonlinearity, Eq. (2) with $\vartheta = 0$, the functions $F$ and $G$ can be calculated explicitly as (see Appendix H)

$$ 
\begin{aligned}
F(h, \sigma^2) &= r_m \Phi(\beta_{\text{eff}}(\sigma^2)h), \\
G(h, \sigma^2) &= r_m^2 \Phi(\beta_{\text{eff}}(\sigma^2)h)
\end{aligned}
\quad (23) $$

$$ -2r_m^2 T \left( \beta_{\text{eff}}(\sigma^2)h, \frac{1}{\sqrt{1 + 2\beta^2\sigma^2}} \right) - F^2(h, \sigma^2) \quad (24) $$

with $T$ being the Owen’s $T$ function and

$$ \beta_{\text{eff}}(\sigma^2) = \frac{\beta}{\sqrt{1 + \beta^2\sigma^2}} $$

being the effective steepness of the sigmoidal nonlinearity at the population level. Because $\beta_{\text{eff}}(\sigma^2) \le \beta_{\text{eff}}(0) = \beta$, the effective nonlinearity that governs the population dynamics has a reduced steepness compared to the single-neuron
<figure><img src="image_6.png" /><figcaption>FIG. 3. Analysis of the deterministic mean-field dynamics. (a) Graphical solution of the fixed-point equation (25), here illustrated for $\mu_0 = 5$ mV and $w = -0.7$ mV s. Solid lines: function $F(\bar{h}, \sigma^2(\bar{h}))$ for MF1 [orange, $\sigma^2(\bar{h}) \equiv 0$] and MF2 [blue, $\sigma^2(\bar{h})$ as in Eq. (25)]. Black dashed line: Left-hand side of Eq. (25). Cyan dash-dotted line: $F$ with the second argument fixed at the fixed-point value $\sigma_0^2$. (b) ($d, w$)-phase diagram with Hopf bifurcation lines (onset of oscillations). Bifurcation lines of MF2 and MFsparselim explain stable asynchronous activity and oscillatory instability. The heat map shows the standard deviation of the population rate $r(t)$ in the quenched network ($p = 0.1$) to indicate the presence of oscillations. Hopf bifurcations, calculated from Eq. (30), are displayed for MF1 (orange line), MF2 (blue line), and MFsparselim (brown line). (c) Sample trajectories of population rate $r(t)$ for quenched network (green), MF1 (orange), MF2 (blue), and MFsparselim (brown) with $d = 1$ ms. [(d)–(f)] Linear response of the time-dependent mean population rate (“rate susceptibility”) to a weak modulation of the external current $\mu(t)$ for quenched network (green circles), MF1 (orange), MF2 (blue), and MFsparselim (dashed brown). The MF theories are given by Eq. (34), where $p$ is set to 1 for MF1, and the sparse limit $N \to \infty$, $pN = C$ = const is taken for MFsparselim. For each panel, top and bottom plots show the amplitude $|\tilde{\phi}_r|$ and the phase $\phi_r$ of the susceptibility, respectively. Dashed horizontal line: $\phi_r = -\pi/2$. (d) $p = 0.1$ ($C = 100$), $d = 0$; (e) $p = 0.4$ ($C = 400$), $d = 1$ ms; (f) $p = 0.95$ ($C = 950$), $d = 0$. Other parameters unless otherwise indicated: $N = 1000$, $w = -1$ mV s, $\beta = 5$ mV$^{-1}$, $\mu_0 = 10$ mV, $\sigma_{\text{ext}} = 0$.</figcaption></figure>

nonlinearity [Fig. 3(a)] and this reduction is caused by the strictly positive spread $\sigma > 0$ of the membrane potentials in a nonfully connected random network.

## IV. ANALYSIS OF THE MESOSCOPIC MODEL

The mesoscopic dynamics MF2, Eq. (20), allows us to calculate and analyze the first- and second-order statistics of mesoscopic quantities such as the stochastic population rate $r(t)$, the population-averaged input potential $\bar{h}(t)$ and the population activity $A_N(t, \Delta t)$. We will mainly focus on the second-order statistics, the analysis of which will be based on a linearization of the dynamics about the fixed points of the noiseless system. Therefore, we will first calculate the fixed points and their stability. Although not the focus of the paper, this analysis will also yield approximations for first-order statistics including the stationary mean population activity and the linear response of the mean population activity to time-dependent external stimuli.

### A. Stationary mean population activity

We first calculate the fixed-points $\bar{h}(t) = \bar{h}_0$, $\sigma^2(t) = \sigma_0^2$, and $\xi(t) = 0$ of the deterministic system when the noise terms in Eq. (20) are set to zero and the external drive is constant, $\mu(t) = \mu_0$. In the following, we assume an inhibitory network ($w < 0$) balanced with a positive external drive $\mu_0 > 0$. In this case, the system exhibits a unique fixed point as shown in the Appendix I. If the fluctuations are sufficiently small and the fixed point is stable, then the quantity $r_0 := F(h_0, \sigma_0^2)$ provides an approximation of the mean stationary population activity (or, equivalently, the mean stationary firing rate of neurons). At equilibrium, the input potential $h_0$ is given by the fixed-point equation

$$ \frac{\bar{h}_0 - \mu_0}{w} = F\left(\bar{h}_0, \frac{w(1-p)}{2\tau pN}(\bar{h}_0 - \mu_0)\right) \quad (25) $$

(see Appendix B). Furthermore, the fixed points of the MF1 are obtained by setting $p = 1$ in Eq. (25) corresponding to setting $\sigma_0^2 = 0$ in the sigmoidal function $F(h_0, \sigma_0^2)$. The solutions for both MF1 and MF2 can be determined graphically [Fig. 3(a)]: For MF1, $\sigma_0^2 = 0$ and therefore the sigmoidal function $F(h_0, 0)$ is steeper compared to MF2, for which $\sigma_0^2$ is strictly positive.

While the fixed points can be obtained numerically, for the analytical theory we use the approximations

$$ r_0 \approx -\frac{\mu_0}{w}, \quad (26) $$

$$ \bar{h}_0 \approx \Phi^{-1}\left(-\frac{\mu_0}{w r_m}\right) \sqrt{\beta^{-2} + \sigma_0^2}, \quad (27) $$

$$ \sigma_0^2 \approx -\frac{w\mu_0(1-p)}{2\tau Np}, \quad (28) $$

for stimuli $0 < \mu_0 < -w r_m$, where the firing rates are low ($0 < r < r_m/2$), see Appendix B. The approximation holds for large $\beta \gg 1/\sigma_0^2$ and for $-\mu_0/w \gg (1-p)/(\tau C)$ and $-\mu_0/w > 0.023 r_m$. The fixed-point solutions yield theoretical predictions for the stationary first-order statistics such as the stationary mean firing rate $\langle r \rangle$, mean membrane potential $\langle \bar{h} \rangle$, and mean membrane-potential spread $\langle \sigma^2 \rangle$. However, we postpone the discussion of these predictions to the Appendix B, because our main focus is the second-order statistics.
## B. Stability of the fixed-point solutions

In the presence of a nonzero transmission delay $d > 0$, the fixed-point solution of the mean-field model may become unstable, leading to an oscillatory regime, see, e.g., Ref. [73]. For such an instability to correctly predict oscillatory behavior of our microscopic network model, our mean-field theory must accurately describe the recurrent fluctuations in the network.

The onset of oscillations can be understood by a linear stability analysis. To this end, we linearize the noiseless system around the fixed point. Small deviations of the mesoscopic variables, $\delta\bar{h} = \bar{h}(t) - \bar{h}_0$, $\delta\sigma^2 = \sigma^2(t) - \sigma_0^2$, obey the linearized dynamics

$$ \frac{d}{dt}X(t) = TX(t) + WX(t-d), \quad (29) $$

where we introduced the deviation vector $X(t) = (\delta\bar{h}(t), \delta\sigma^2(t), \xi(t))^T$ and the matrices

$$ T = -\frac{\operatorname{diag}(1, 2, 1)}{\tau}, \quad W = \left(\frac{w}{\tau}, \frac{w^2(1-p)}{\tau^2 p N}, 0\right)^T L^T $$

with $L^T = (F_h, F_\sigma, N^{-1/2})$. Here $F_h$ and $F_\sigma$ denote the partial derivatives of $F(h, \sigma^2)$ with respect to $h$ and $\sigma^2$ at the fixed point (exact and approximate expressions for $F_h$ and $F_\sigma$ are given in the Appendix D). We note that any perturbation in $\xi$ exponentially decays to zero with a time constant $\tau$. It is therefore sufficient to look for solutions of the form $X(t) = e^{\lambda t}(\hat{h}(\lambda), \hat{\sigma}(\lambda), 0)^T$ with some constant parameter $\lambda$. The sign of the real part of $\lambda$ then determines the stability of the fixed point; in particular if $\operatorname{Re}(\lambda) > 0$ the fixed point is unstable. As shown in Appendix E, the parameter $\lambda$ must satisfy the characteristic equation

$$ \lambda\tau = -1 + wF_h e^{-\lambda d} + \frac{wF_h \hat{F}_\sigma e^{-2\lambda d}}{2 + \lambda\tau - \hat{F}_\sigma e^{-\lambda d}}, \quad (30) $$

where $\hat{F}_\sigma := w^2 \frac{1-p}{\tau C} F_\sigma$.

To find the onset of oscillations, we look for a Hopf bifurcation, at which $\lambda = i\omega$ with a frequency $\omega \in \mathbb{R} \setminus \{0\}$. Imposing such purely imaginary values of $\lambda$ in Eq. (30) allows us to find the “Hopf boundary” of the oscillatory phase in the parameter space. In the following, we investigate the stability depending on the transmission delay $d$ and the coupling strength $w$ [Fig. 3(b)]. In the $(d, w)$-space, we can see that the bifurcation lines of the noiseless systems significantly differ between second-order and first-order MF dynamics. For sufficiently strong coupling strength, the first-order MF dynamics exhibits a Hopf boundary at small values of $d$, predicting oscillations already at small delays. In contrast, the second-order theory predicts an oscillation onset for much larger delays. Interestingly, MF2 predicts that a stronger coupling strength $w$ requires a larger delay to enable oscillations, whereas the first-order model requires a smaller delay. We note that the bifurcation line obtained for MFsparselim is close to MF2. The bifurcation lines qualitatively explain the presence or absence of large noisy oscillations in the mean-connectivity and random network models [Figs. 3(b) and 3(c)]. Furthermore, the simulations of the MF dynamics indicate that the Hopf bifurcation is supercritical: Varying the delay, the variance of the population rate $r$ changes continuously at the bifurcation and no bistability between fixed point and limit cycle is observed.

In conclusion, our stability analysis of MF2 shows that connectivity disorder leads to a stabilization of the nonoscillatory stationary state. In particular, the disordered connectivity leads to incoherent fluctuations of the recurrent input, which flattens the effective nonlinearity of the macroscopic dynamics and thus changes the stability properties of the fixed point.

## C. Linear response of the mean activity

We now analyze how the network responds to a small, time-dependent perturbation of the external current, $\mu(t) = \mu_0 + \mu_1(t)$ by calculating the linear response function [74]. To this end, we linearize our mesoscopic mean-field dynamics, Eq. (20), around the fixed point and obtain a generalization of Eq. (29):

$$ \frac{d}{dt}X(t) = TX(t) + WX(t-d) + M(t) + B\xi(t). \quad (31) $$

Here the vector $X(t)$ represents again the deviation from the fixed point as in the previous section and the three-dimensional vector $M(t)$ represents a general small perturbation to the system. Furthermore, the last term represents the finite-size noise, where

$$ B = \begin{pmatrix} \frac{w}{\tau} \sqrt{\frac{r_0}{N}} & 0 \\ 0 & 0 \\ 0 & \sqrt{\frac{2G_0}{\tau}} \end{pmatrix}, \quad \xi(t) = \begin{pmatrix} \eta(t) \\ \zeta(t) \end{pmatrix} $$

and $G_0 = G(\bar{h}_0, \sigma_0^2)$. In the linearized dynamics we neglected terms of order $O(1/N)$ and higher. In the Fourier domain, the solution reads

$$ \bar{X}(\omega) = \tilde{\chi}(\omega)(\tilde{M}(\omega) + B\tilde{\xi}(\omega)), \quad (32) $$

where we introduced the susceptibility matrix

$$ \tilde{\chi}(\omega) = [i\omega I_3 - T - W e^{-i\omega d}]^{-1}. \quad (33) $$

Here $I_3$ denotes the 3×3 identity matrix and the tilde notation $\tilde{f}(\omega)$ for a function $f(t)$ denotes its Fourier transform, i.e., $\tilde{f}(\omega) = \int_{-\infty}^{\infty} f(t)e^{-i\omega t} dt$. The explicit expressions for the elements of $\tilde{\chi}$ are given in Appendix F, Eq. (F1).

The susceptibility matrix provides the response of the ensemble mean $\langle X(t) \rangle$ of the deviation from the fixed point of the full nonlinear system to an infinitesimally small perturbation $M(t)$ in the Fourier domain via the relation $\langle \tilde{X}(\omega) \rangle = \tilde{\chi}(\omega) \tilde{M}(\omega)$. Similarly, the response of the ensemble-averaged firing rate $\langle r(t) \rangle = F(h_0, \sigma_0^2) + \delta r(t)$ is given in the Fourier domain by $\delta r(\omega) = L^T \tilde{\chi}(\omega) \tilde{M}(\omega)$. For an external stimulus $\mu(t) = \mu_0 + \mu_1(t)$, the perturbation is $\tilde{M}(\omega) = (\tilde{\mu}_1(\omega)/\tau, 0, 0)^T$. In this case, we then have $\delta r(\omega) = \tilde{\chi}_r(\omega) \tilde{\mu}_1(\omega)$, where we introduced the rate susceptibility

$$ \tilde{\chi}_r(\omega) = \frac{1}{\tau} [F_h \tilde{\chi}_{11}(\omega) + F_\sigma \tilde{\chi}_{21}(\omega)]. \quad (34) $$

Here we used that $\tilde{\chi}_{31} \equiv 0$.

The rate susceptibility is in general complex-valued with absolute value $|\tilde{\chi}_r(\omega)|$ and phase $\phi_r(\omega)$. A straightforward interpretation is that the system driven by a sinusoidal external input modulation $\mu_1(t) = \epsilon \sin(\omega t)$ with small amplitude $\epsilon$ will respond with a sinusoidal rate modulation with amplitude
$\epsilon|\tilde{\chi}_r(\omega)|$ and a phase shift $\phi_r(\omega)$. For all models studied here, the amplitude decreases for high frequencies with a power law with exponent $-1$ [Figs. 3(d)-3(f)]. The MF2 theory agrees very well with the quenched microscopic network, in contrast to the MF1 theory which massively overestimates the amplitude at high frequencies. The quenched network shows a slightly smaller amplitude than predicted by the second-order MF theory but is still reasonably well approximated by the latter. While the amplitude vanishes at low frequencies as expected, the phase shift approaches $-\frac{\pi}{2}$ in the limit of large frequencies [Figs. 3(d)-3(f)]. We can understand this behavior analytically by considering the limits of Eq. (34). For $\omega = 0$ we obtain

$$ \tilde{\chi}_r(0) = \frac{F_h}{1 - wF_h - \frac{w^2}{2\tau}F_\sigma\left(\frac{1}{C} - \frac{1}{N}\right)}, \quad (35) $$

which is real valued and has therefore a phase zero. The amplitude in the low-frequency limit is approximately the negative inverse of the coupling strength $-w^{-1}$, under the condition that $1 \ll wF_h$ and $C$ large. Similarly, for large frequencies we find asymptotically

$$ \tilde{\chi}_r(\omega) \sim -i \frac{F_h}{\tau} \cdot \frac{1}{\omega}, \quad \omega \to \infty. \quad (36) $$

In this limit, the phase converges to $-\frac{\pi}{2}$ and the amplitude follows a power law $\omega^{-1}$ scaled by $F_h/\tau$. The different values for the derivative $F_h$ explain the difference between the MF1 and the MF2 theory. To understand this difference analytically, we use the following approximation for the slope of the transfer function at the fixed point (see Appendix D):

$$ F_h \approx \frac{r_m Q\left(-\frac{\mu_0}{w r_m}\right)}{\sqrt{2\pi(\beta^{-2} + \sigma_0^2)}}, \quad (37) $$

where $\sigma_0^2$ is given by Eq. (28). Furthermore, we introduced the nondimensional function

$$ Q(x) := \exp\left[-\frac{1}{2}\Phi^{-1}(x)^2\right] \quad (38) $$

defined on the interval $[0, 1]$. This function has an inverted-U shape and is symmetric with respect to $x = \frac{1}{2}$ [a graph of $Q$ is shown in Fig. 6(c), dotted lines]. Equation (37) clearly shows that the factor $F_h$ increases with the connection probability $p$ because $\sigma_0^2$ decreases with $p$, cf. Eq. (28). The MF1 theory corresponding to the mean-connectivity network (without additional noise) is obtained by setting $p = 1$ in Eq. (28). In this case, the dispersion $\sigma_0^2$ vanishes, and hence the slope $F_h$ is maximized. Therefore, the MF1 and MF2 theories predict the same power-law behavior of the rate susceptibility with $\omega$ at large frequencies, but with a higher magnitude for the first-order model [Fig. 3(d)]. With increasing connection probability $p$ the susceptibility approaches the first-order MF theory. The low-frequency limit is in both mesoscopic models almost independent of $p$ as the low-frequency limit is roughly given by $-w^{-1}$. Simulations of the annealed network network (not shown in the figure) match the predictions of the MF2 theory while the quenched network model has slightly lower amplitude and slightly lower phase for intermediate frequencies. Overall, the MF2 theory yields a strongly improved prediction compared to the MF1 theory.

The synaptic delay $d$ appears in $\tilde{\chi}_r$ only via the complex exponential in Eq. (33). On the other hand, the delay does

not influence the low-frequency limit and the high-frequency asymptotics, Eqs. (35) and (36) of the response, respectively. For intermediate frequencies the complex exponential causes a modulation and a peak in the amplitude as a function of $\omega$ [Fig. 3(e)]. Simulations of the annealed network model perfectly confirm (not shown in figure) the prediction of the second-order MF theory, whereas the quenched network model has a less pronounced peak at a lower frequency. The susceptibility of MFsparselim fits well with the quenched network for small or medium connection probability, $p = 0.1$ and $p = 0.4$ [Figs. 3(d) and 3(e)]. For larger $p$, however, there is a noticeable discrepancy for MFsparselim, while MF1 and MF2 have a good match with the quenched network [at $p = 0.95$, see Fig. 3(f)]. In marked contrast to MF1 and MFsparselim, the predictions of MF2 are good over the entire range of connection probabilities $p$.

### D. Power spectral density

Now we turn to the second-order statistics of the stationary dynamics. We assume that the external stimulus $\mu(t)$ is a stationary process with given mean $\mu_0$ and power spectral density $S_{\mu\mu}(\omega)$. We begin with the power spectrum of the population activity $A(t)$ and population rate $r(t)$. This statistics has been calculated analytically using linear-response theory for finite-size networks with uniform full connectivity [47,71]. In a similar fashion, we here derive the power spectrum via a linearization of our second-order MF dynamics. To this end, we start with the power spectral density matrix $S(\omega)$ of the mean-field state vector $X(t)$ [cf. Eq. (29)] defined via the relation

$$ \langle \tilde{X}(\omega) \tilde{X}^*(\omega') \rangle = 2\pi S(\omega)\delta(\omega - \omega'), \quad (39) $$

where * denotes the conjugate transpose of a matrix. Using Eq. (32) for the linearized system, we find the following approximation for $S(\omega)$ describing weak fluctuations around the fixed points (see Appendix G):

$$ S_{ij}(\omega) = \frac{1}{\tau^2} \left[ S_{\mu\mu}(\omega) + \frac{w^2 r_0}{N} \right] \tilde{\chi}_{i1}(\omega) \tilde{\chi}_{j1}^*(\omega) + \frac{2G(\tilde{h}_0, \sigma_{h,0}^2)}{\tau} \tilde{\chi}_{i3}(\omega) \tilde{\chi}_{j3}^*(\omega). $$

The matrix $S$ gives us the (cross-) spectral densities among $\tilde{h}$, $\sigma^2$, and $\tilde{\xi}$. For the power spectral densities of the population activity $A(t)$, we use the linear approximation

$$ \tilde{A} \approx \tilde{r} + \sqrt{\frac{r_0}{N}} \tilde{\eta}, \quad \tilde{r} \approx F_h \tilde{h} + F_\sigma \tilde{\sigma}^2 + \frac{1}{\sqrt{N}} \tilde{\xi}, \quad (40) $$

which yields approximate expressions for the power spectrum of the stochastic population rate $r(t)$,

$$ S_{rr}(\omega) \approx F_h^2 S_{11}(\omega) + F_\sigma^2 S_{22}(\omega) + \frac{1}{N} S_{33}(\omega) \\ + 2F_h F_\sigma \mathfrak{Re}(S_{12}(\omega)) \\ + 2 \frac{F_h \mathfrak{Re}(S_{13}(\omega)) + F_\sigma \mathfrak{Re}(S_{23}(\omega))}{\sqrt{N}}, \quad (41) $$
<figure><img src="image_7.png" /><figcaption>FIG. 4. Power spectral densities of the population activity and population rate. (a) Power spectrum of $A_N(t; \Delta t)$ for quenched (green circles) and annealed network (red squares) and first- and second-order theory [orange and blue lines, Eqs. (G1) and (42), respectively]. (b) Same for the power spectrum of $r(t)$ [theory: Eq. (41)] Parameters: $w = -1$ mV s, $N = 500$, $p = 0.07$, $\mu(t) = \mu_0 = 10$ mV, $\sigma_{\text{ext}} = 0$, $\Delta t = 0.1$ ms.</figcaption></figure>

and of the population activity $A(t)$:

$$S_{AA}(\omega) \approx S_{rr}(\omega) + \frac{r_0}{N} + 2\sqrt{\frac{r_0}{N}}\Re(S_{r\eta}(\omega)). \quad (42)$$

The cross spectral density $S_{r\eta}$ that appears here can be calculated using Eqs. (40) and (32):

$$S_{r\eta}(\omega) \approx \frac{w}{\tau} \sqrt{\frac{r_0}{N}} [F_h \tilde{\chi}_{11}(\omega) + F_\sigma \tilde{\chi}_{21}(\omega)].$$

In Fig. 4, we compare the theoretical predictions, Eqs. (41) and (42), to the corresponding power spectra obtained from simulations of the quenched and annealed network. In all cases, the power spectra of the population activity display a pronounced trough at low frequencies [Fig. 4(a)]. The MF1 theory severely underestimates the observed low-frequency power, whereas the MF2 theory matches well the microscopic model under the annealed approximation. Nonetheless, even with this marked improvement over MF1, the annealed network's low-frequency power remains lower than that of the quenched network, exposing a limitation of the annealed approximation. The discrepancy of the power spectra of quenched and annealed networks is caused by the partial elimination of the effect of the single-neuron autocorrelations on the recurrent input when applying the annealing (i.e., the permanent resampling of connectivity). We note that the elimination of temporal correlations only concerns the incoherent part of the recurrent input, as expressed by the white noise in Eq. (12). In contrast, for finite $N$, the coherent part of the recurrent input fluctuations (finite-size fluctuations of population activity) still contains effects of spike-train autocorrelations, even in the annealed network. Finally, we mention that we have not included MFsparselim in Fig. 4 because its stationary dynamics is characterized by a

<figure><img src="image_8.png" /><figcaption>FIG. 5. Stationary variance of the population rate and mean input potential for varying connection probability $p$. (a) The number of neurons is fixed to $N = 1000$ (hence $C = 1000p$). (b) The in-degree is fixed to $C = 100$ (hence $N = 100/p$). Left panels of (a) and (b): Variance of the population rate $r(t)$ [Symbols: simulations of microscopic networks as indicated in the legend; orange and blue solid lines: full theory Eq. (43) and Eqs. (G4); black dashed line: analytical approximation for the variances based on the second-order MF theory, Eqs. (46), (47) and Eqs. (45), (49) for $\text{var}(r)$ and $\text{var}(\bar{h})$, respectively; dash-dotted line: MFsparselim.] Right panels of (a) and (b): Same as left panels but for variance of $\bar{h}$. Parameters: $\mu_0 = 10$ mV, $w = -1$ mV s, $\sigma_{\text{ext}} = 0$ mV, $\beta = 5$ mV$^{-1}$.</figcaption></figure>

deterministic equilibrium without any fluctuations and a vanishing power spectrum.

### E. Variances

The stationary variances of $r(t)$ and $\bar{h}(t)$ can be calculated by standard methods for the linearized system in the absence of synaptic delay, $d = 0$ (see Appendix G). Furthermore, we assume that the external stimulus $\mu(t)$ has mean $\mu_0$ and white Gaussian fluctuations of strength $\sigma_{\text{ext}}^2$ as in Sec. II B. Under these assumptions, one eventually obtains

$$\text{var}(r) \approx F_h^2 \sigma_{11} + F_\sigma^2 \sigma_{22} + 2F_h F_\sigma \sigma_{12} + \frac{2F_h}{\sqrt{N}} \sigma_{13} + \frac{2F_\sigma}{\sqrt{N}} \sigma_{23} + \frac{1}{N} G(\bar{h}_0, \sigma_0^2). \quad (43)$$

where $\sigma_{ij} = \langle X_i(t)X_j(t) \rangle$ are the elements of the covariance matrix of the state vector $X(t)$. In particular, $\sigma_{11}$ yields the variance of the population-averaged input potential $\bar{h}(t)$. Approximate analytical expressions for $\sigma_{ij}$ can be calculated using the linearized system, Eq. (31), and are given in the Appendix G. Our theory, Eq. (43), agrees well with microscopic simulations (Figs. 5 and 6), discussed in more detail below. In general, we find a good quantitative match between the theory and simulations of the annealed network. While, in the case of quenched random connectivity, the theory again exhibits quantitative deviations, our second-order MF theory still captures the parameter dependence of the variance qualitatively and yields a much better quantitative prediction than the first-order MF theory. Again, MFsparselim without any
<figure><img src="image_9.png" /><figcaption>FIG. 6. Variance of the stochastic population rate $r(t)$ for varying mean external stimulus $\mu_0$. (a) Moderate-size network of $N = 10^3$ neurons whose inputs exhibit internally generated, finite-size noise but no external noise, $\sigma_{\text{ext}} = 0$. Different in-degree $C$ are encoded by color. Left: Response to a step stimulus, where $\mu_0$ jumps from 28 to 50 mV at time $t = 0.04$ s. Blue and orange lines: simulations of MF2 for $C = 100$ ($p = 0.1$) and $C = 1000$ ($p = 1$), the latter ($p = 1$) equivalent to MF1 or the mean-connectivity approximation of the former ($p = 0.1$). Black dotted: Corresponding theory for the stationary variances, Eq. (43). Gray dashed: Stimulus onset. Inset: Magnified data for $C = 100$. Symbols represent stationary variances of network simulations and theory, see description of right panel. Right: Stationary variance of $r(t)$ vs mean stimulus strength $\mu_0$. Symbols: Simulations for $C = 100$ of quenched network (green circles), MF1 (orange triangles), and MF2 (blue diamonds). Solid lines: Full MF1 (orange) and MF2 theory (blue). Black dashed: Corresponding analytical approximations, Eq. (47); dash-dotted: MFsparselim. The gray dashed line marks the border of the region of interest $\mu_0 \le -wr_m/2$ where the hazard function at the fixed point, $\phi(\bar{h}_0)$ is convex. (b) Analogous to (a) for a much larger network ($N = 5 \times 10^4$) with negligible internally generated finite-size noise but with additional external noise $\sigma_{\text{ext}}^2 = 1$ mV$^2$. The connection probability is varied by different values of $C$ (encoded by color). Orange line: Common MF1 theory or mean-connectivity network. Dash-dotted lines: MFsparselim. (c) Slope $\hat{F}_h$ vs mean stimulus $\hat{\mu}_0$. Dotted: $p = 1$ corresponding to the function $Q$; solid line: $p = 0.7$. Gray dashed as in (a) and (b) right. (d) Same for the product $\hat{\mu}_0 \hat{F}_h$. [(c) and (d)] $N = 1000$, $\hat{\beta} = 500$ [equivalent to $\beta = 5$ mV$^{-1}$ in (a) and (b)].</figcaption></figure>

noise in the external input deterministically converges into an equilibrium with zero variance.

1. Dependence on the connection probability

Which connection probability causes the largest variability in our network model? To address this question, we analyze the variances of the population rate $r$ and the average input potential $\bar{h}$ as a function of $p$ (Fig. 5). To avoid a trivial dependence through the mean firing rate, we vary $p$ such that the mean firing rates remain roughly unchanged. According to Eq. (26), a constant rate can be achieved by keeping the coupling strength $w = NpJ$ constant. In an experimental setting, where the synaptic efficacy $J$ is the basic physiological parameter, this constraint could be realized in two different scenarios, either by fixing the network size $N$ and changing the synaptic efficacy $J$ such that $Jp = \text{const.}$ or by fixing both the number of input connections $C = Np$ and the synaptic efficacy $J$.

In the first scenario with constant $N$, we find that the variance of the rate $r(t)$ increases with connection probability in microscopic network simulations [Fig. 5(a), left], although the variance of the input potential $\bar{h}$ decreases [Fig. 5(a), right]. MF1 cannot explain this behavior because the model equation (21) does not depend on the parameter $p$ when $w$ is fixed and $N$ is constant. Thus, for all $p \in [0, 1]$, the variance predicted by the first-order MF theory corresponds to the fully connected case $p = 1$, i.e., to the mean-connectivity network. In contrast, the dependence on $p$ of the variances measured from simulations of the annealed network are quantitatively well reproduced by the second-order MF theory, Eq. (43). The second-order theory also captures the variances of the quenched network qualitatively. The quantitative agreement is also reasonably good; however, we find that the rate variance is overestimated and the variance of the population-averaged input potentials is underestimated by the annealed network. Nevertheless, the prediction of variances by the second-order theory yields a significant improvement over the first-order MF theory.

In order to gain an analytical understanding, we use the following explicit approximations of the variances (see Appendix G):

$$ \text{var}(\bar{h}) \equiv \sigma_{11} \approx \frac{-w\mu_0 + \sigma_{\text{ext}}^2}{2\tau N}, \quad \text{var}(r) \approx F_h^2 \text{var}(\bar{h}). \quad (44) $$

In the last approximation, we have used Eq. (43) and the fact that, for the parameters considered in this study, the variance of $r$ is vastly dominated by the first term in that equation. By contrast, the other terms in Eq. (43) are small enough to be ignored. Equation (44) reveals that an increasing slope $F_h$ decreases the variance of $\bar{h}$ in an inhibitory network ($w < 0$). Intuitively, this effect can be understood by the increase of the effective leakage $(1 - wF_h)/\tau$ in the linearized dynamics of $\bar{h}$ which dampens the input fluctuations. In contrast, the variance of the population rate can increase with increasing slope because the quadratic prefactor $F_h^2$ can overcompensate the decrease of $\text{var}(\bar{h})$. To make this argument more clear in the theoretical discussion, we further assume that $-wF_h \gg 1$ (which holds for our parameter setting). Under this assumption, the variance of the population-averaged input potentials
becomes inversely proportional to $F_h$,

$$ \mathrm{var}(\bar{h}) \approx \left( \frac{\mu_0}{2\tau N} + \frac{\sigma_{\mathrm{ext}}^2}{-2w} \right) \frac{1}{F_h}. \quad (45) $$

In turn, the quadratic factor $F_h^2$ due to the conversion from variance of $\bar{h}$ to variance of $r$ leads effectively to a linear dependence of $\mathrm{var}(r)$ on $F_h$,

$$ \mathrm{var}(r) \approx \left( \frac{\mu_0}{2\tau N} + \frac{\sigma_{\mathrm{ext}}^2}{-2w} \right) F_h. \quad (46) $$

The variance formulas (45) and (46) share a common factor that represents the sum of two contributions to the variability. The first term of the sum corresponds to the intrinsically generated variability caused by finite-size fluctuations of order $\sqrt{r_0/N}$ in the mesoscopic MF dynamics, Eq. (20). We recall that the dependence of these fluctuations on the firing rate $r_0 = -\mu_0/w$, arises from the Poisson spiking noise at the microscopic level. Thus, the Poisson property that the variance scales proportionally to the mean causes finite-size fluctuations that increase with the mean stimulus strength $\mu_0$ (see below). We emphasize that this contribution to the variability is a clear finite-size effect that vanishes in the limit $N \to \infty$. In contrast, the second term of the sum in Eqs. (45) and (46) represents externally generated variability proportional to the variance of the common external stimulus.

For constant $N$, the $p$ dependence in both formulas, Eqs. (45) and (46), is solely contained in the slope $F_h$ of the transfer function. As discussed in Eq. (37) above, the slope $F_h$ is a monotonously increasing function of $p$. The slope dependence thus explains the increase of the variance of the population rate [Fig. 5(a), left] and the decrease of the variance of the mean input potential [Fig. 5(a), right], respectively.

In the second scenario with constant $C$ and varying population size $N = C/p$, we observe different dependencies of the variances on $p$ compared to the first scenario with constant $N$ [Fig. 5(b)]. In simulations of microscopic networks and second-order MF dynamics, the rate variance increases supra-linearly with $p$ and the variance of the mean input potential $\bar{h}$ exhibits a nonmonotonic behavior with a maximum inside the range $0 < p < 1$. Again, the first-order MF theory does not capture these behaviors. Note that there is a good qualitative agreement between second-order MF theory and microscopic network simulations. As in the first scenario, however, the variance of the average input potential in the quenched network is underestimated by the annealed network, which, in turn, is quantitatively well matched by the second-order MF theory.

To gain a theoretical understanding of these observations, we use our analytical approximation Eq. (46). For constant $C$, the approximation can be rewritten as

$$ \mathrm{var}(r) \approx \frac{r_m^2 \hat{\beta}}{\sqrt{8\pi}} \left( \frac{\hat{\mu}_0 p}{\hat{\tau} C} + \hat{\sigma}_{\mathrm{ext}}^2 \right) \hat{F}_h, \quad (47) $$

where we have introduced the dimensionless parameters $\hat{\mu}_0 = -\mu_0/(wr_m)$, $\hat{\sigma}_{\mathrm{ext}}^2 = \sigma_{\mathrm{ext}}^2/(wr_m)^2$, $\hat{\beta} = -\beta wr_m$, $\hat{\tau} = \tau r_m$, as well as the dimensionless slope of the transfer function

$$ \hat{F}_h := \sqrt{2\pi} \frac{F_h}{\beta r_m} \approx \frac{Q(\hat{\mu}_0)}{\sqrt{1 + \frac{\hat{\beta}^2(1-p)}{2\hat{\tau}C} \hat{\mu}_0}} \quad (48) $$

[cf. Eq. (37)]. The expression for the rate variance in Eq. (47) is a product of two monotonously increasing functions of $p$—a linear function (the prefactor of $\hat{F}_h$) and a strictly convex function ($\hat{F}_h$). Hence, the variance of the population rate is a monotonously increasing, strictly convex function of $p$ explaining our observation in Fig. 5(b) (left). In contrast, the variance of the mean input potential, Eq. (45), can be rewritten for constant $C$ and dimensionless parameters as

$$ \mathrm{var}(\bar{h}) \approx \frac{\sqrt{2\pi} w^2 r_m^2}{2\hat{\beta}} \left( \frac{\hat{\mu}_0 p}{\hat{\tau} C} + \hat{\sigma}_{\mathrm{ext}}^2 \right) \frac{1}{\hat{F}_h}. \quad (49) $$

This expression is a product of a monotonously increasing and a monotonously decreasing function of $p$ (the prefactor of $1/\hat{F}_h$ and $1/\hat{F}_h$, respectively). The product may therefore exhibit a nonmonotonic behavior. A closer inspection of Eq. (49) reveals a maximum at an intermediate connection probability

$$ p_{\max} = \frac{2}{3} \left( 1 + \frac{2\tau C}{-w\mu_0} \left( \beta^{-2} - \frac{\sigma_{\mathrm{ext}}^2}{4} \right) \right) $$

in the interval $(0, 1)$ if $-2/3 < \frac{2\tau C}{-w\mu_0}(\beta^{-2} - \sigma_{\mathrm{ext}}^2/4) < 1/2$. In the limit of a steep single-neuron transfer function ($\beta \to \infty$) and in the absence of external noise, the maximum is attained at $p_{\max} = 2/3$. Such maximum is indeed verified by simulations [Fig. 5(b), right]. Furthermore, the variance $\mathrm{var}(\bar{h})$ approaches zero for $p \to 0$ and approaches a nonvanishing value for $p \to 1$.

In contrast to the second-order theory, the first-order MF theory cannot explain the strictly convex and nonmonotonic behavior of the variances of the population rates and input potentials, respectively [Fig. 5(b)]. In the second scenario with constant $C$, the first-order MF theory is obtained from the above formulas by fixing the ratio $p/C$ and letting $p \to 1$. Thus, the common factor $(\frac{\hat{\mu}_0 p}{\hat{\tau} C} + \hat{\sigma}_{\mathrm{ext}}^2)$ in Eqs. (47) and (49) remains unchanged, whereas the second factor becomes independent of $p$. Because the common factor is an increasing linear function of $p$, the variances are thus also increasing linear functions, in line with simulations of the first-order mean-field model (but in contradiction to microscopic simulations and second-order MF theory).

## 2. Dependence on the mean stimulus strength

We started our paper with network simulations pointing out marked quantitative and qualitative discrepancies of the mean-connectivity network with respect to the stimulus dependence of the rate variability (Fig. 1). In particular, we observed a strong overestimation of the rate variance by the mean-connectivity network compared to the original quenched network. At first sight, it may appear counterintuitive that more homogeneous connectivity leads to a stronger variance in the population rate. This apparent paradox results from the fact that in mean-connectivity networks, finite-size fluctuations are “directly” picked up by the steep single-unit transfer function that also determines the population rate, $F(h, \sigma^2) = \phi(h)$ for $\sigma^2 = 0$, whereas the effect of these fluctuations at the population level is attenuated for networks with heterogeneous connectivity due to the flattening of the effective population rate transfer function; the difference in the population rate variance can then be seen from Eq. (46).
We will make this argument more explicit below. We also note that the large variance of the population rate in the mean-connectivity network (or MF1) can be attributed to an increased high-frequency power [Fig. 4(b)].

Furthermore, for a noisy external stimulus and large network size ($N = 50000$), we observed a suppression of the rate variability by an increase in the mean stimulus strength in the quenched network. In contrast, the mean-connectivity network incorrectly predicted an increase in rate variability rather than a suppression. As we will see now, these observations can also be well explained by our second-order MF theory.

Let us study two opposing scenarios: First, we consider a small ("mesoscopic") network without noise in the external stimulus ($\sigma_{\text{ext}} = 0$). In this case, the rate variability solely originates intrinsically from the finite-size fluctuations. Second, we consider a large ("macroscopic") network with a noisy external stimulus ($\sigma_{\text{ext}} > 0$). In this case, the rate variability is purely externally generated. In both cases, we find that the first-order MF theory largely overestimates the variance of the population rate in the quenched and annealed network, whereas the second-order MF theory, Eq. (43), correctly predicts the magnitude of the rate variance [Figs. 6(a) and 6(b)]. Because the first-order MF dynamics is the exact MF dynamics for the mean-connectivity network, this large deviation explains our observation in Fig. 1.

Furthermore, in the two different scenarios, we find a different monotonicity of the rate variance if we restrict ourselves to the biologically interesting range where the fixed point is in the convex part of the sigmoidal hazard function, $\tilde{h}_0 < 0$ and $r_0 < r_m/2$, i.e., below the inflection point of $\phi(h)$ [cf. Fig. 3(a)]. The upper boundary of this range given by the inflection point corresponds to a single-neuron firing rate at half maximum, $r_m/2$, and hence an external stimulus $\mu_0^* = -wr_m/2$ [Figs. 6(a) and 6(b) right, dashed vertical line]. In the first scenario without external noise, the finite-size induced rate variability increases monotonically for $\mu_0 < \mu_0^*$ for both first- and second-order MF theory. In contrast, in the second scenario with external noise and without finite-size fluctuations, the variance of the rate shows a nonmonotonic behavior with a maximum below $\mu_0^*$ if $p < 1$.

These observations can be understood analytically using the approximation Eqs. (47) and (48). Let us first discuss how the factor $\hat{F}_h$ in Eq. (47) depends on $\hat{\mu}_0$. For $p = 1$, i.e., in the first-order MF theory, $\hat{F}_h$ is exactly given by the function $Q(\hat{\mu}_0)$ displayed by the dotted line in Fig. 6(c). This function, and thus $\hat{F}_h$, grows monotonically for $\hat{\mu}_0 > 0$ until it reaches a maximum at $\hat{\mu}_0^* = -\mu_0^*/(wr_m) = 0.5$, which corresponds to the upper bound $r_m/2$ of the biologically relevant dynamical range. For $p < 1$, the factor $\hat{F}_h$ is the quotient of the symmetric function $Q$ and an increasing function, which shifts the maximum of $\hat{F}_h$ to a value $< 1/2$ [Fig. 6(c), solid line]. Hence, $\hat{F}_h$ exhibits a nonmonotonic behavior in the biologically relevant range $0 < \mu_0 < \mu_0^*$. Besides the shift of the maximum to smaller values of the mean stimulus strength $\hat{\mu}_0$, the maximum also decreases when the connection probability is lowered from $p = 1$ to smaller values. This decrease is a consequence of dividing the function $Q(\hat{\mu}_0)$ by the square-root term that is strictly larger than one for $p < 1$. Importantly, $\hat{F}_h$ can decrease by a large factor if the single-neuron

transfer function is steep, i.e., $\hat{\beta} \gg 1$, which explains the drastic decrease of the rate variance from first- to second-order MF theory [Figs. 6(a) and 6(b), where $\hat{\beta} = 500$, corresponding to $\beta = 5 \text{ mV}^{-1}$ as in Fig. 1].

The behavior of $\hat{F}_h$ just described explains the variance of the population rate for an infinitely large network with external noise [second scenario, Fig. 6(b)]. In fact, the prefactor of $\hat{F}_h$ in Eq. (47) becomes independent of $\hat{\mu}_0$ in the limit when $N = C/p \to \infty$, and hence the dependence on $\hat{\mu}_0$ is fully captured by $\hat{F}_h$. In particular, the stark difference in magnitude of var($r$) between the first- and second-order MF prediction, as well as the qualitative response of the rate variability (suppression or amplification) to increasing stimuli, is explained by how the monotonicity and magnitude of the slope of the population transfer function $\hat{F}_h$ varies with $p$.

In contrast, in the first scenario when $\sigma_{\text{ext}} = 0$ and $N = C/p$ is finite, Eq. (47) tells us that the variance of the population rate is proportional to

$$ \hat{\mu}_0 \hat{F}_h = \frac{\hat{\mu}_0}{\sqrt{1 + \frac{\hat{\beta}^2(1-p)}{2\hat{\tau}C} \hat{\mu}_0}} Q(\hat{\mu}_0). $$

The right-hand side is a product of a monotonically increasing function of $\hat{\mu}_0$ and the function $Q(\hat{\mu}_0)$, which is maximized at $\hat{\mu}_0^* = 1/2$. Therefore, the maximum of $\hat{\mu}_0 \hat{F}_h$ is shifted to values $\hat{\mu}_0 > 1/2$ for all $p \in [0, 1]$ [Fig. 6(d)]. As a consequence, the variance of the population rate var($r$) is a monotonically increasing function of the mean stimulus strength in the biologically relevant range $0 \le \hat{\mu}_0 \le 1/2$ [Fig. 6(a), right $\mu_0 < 50 \text{ mV}$], and thus an increase of the mean stimulus strength always increases the rate variability in the absence of external common noise [Fig. 6(a), left].

To keep the arguments simple, we studied here two clear-cut cases—one with and another without external noise. According to our theory, Eq. (46), the rate variance is a sum of two terms, where each term on its own corresponds to one of these cases. In general, the theory also applies to the mixed case, where the rate variance is the sum of externally and intrinsically induced variability. For example, in the simulations of Fig. 1, we also compared two networks of size $N = 1000$ and $N = 50000$ but the external common noise was present in both cases. Nonetheless, the mesoscopic case of $N = 1000$ with external noise can still be qualitatively understood by the limit case $\sigma_{\text{ext}} = 0$ because the finite-size fluctuations dominate over the external noise. Last, we remark that as $\hat{\mu}_0 \to 0$, e.g., for strongly inhibition-dominated networks, the population rate variance decreases for both the mean-connectivity network (MF1) and MF2 in both cases as can be seen from Figs. 6(c) and 6(d).

## V. CONCLUSIONS

In this work, we set out to derive a stochastic population dynamics for finite-size networks of spiking neurons that correctly accounts for fluctuations at the mesoscopic scale when the connectivity is random—as in the case of nonfull connectivity ($p < 1$). While such regime is biologically highly relevant, to our knowledge it has not yet been systematically investigated within the frameworks of previously developed mean-field theories [29,30]. The low-dimensional, stochastic
population model derived here allows one to study specifically
how the combination of quenched disorder and stochastic
neuronal firing affects the mesoscopic activity of finite-size
neuronal populations. We thus hope to close an important gap
in the use cases of simplified descriptions of networks and to
provide the foundations for a mean-field modeling framework
of neural variability at the mesoscopic population level.

Based on a Poisson assumption for neural firing and a
"dynamically annealed" description of the connectivity, we
approximated the recurrent input to each neuron by its (pop-
ulation) mean, a coherent fluctuating part shared among
all neurons, and an individually fluctuating contribution re-
flecting the connectivity disorder. Eventually, we obtained
three coupled stochastic differential equations (SDEs) that
describe the evolution of the mean input potential $\bar{h}$ within
the population, the input potential's variance $\sigma^2$, as well as
a colored noise $\xi$ in the population rate that vanishes in
both limits $N \to \infty$ and $p \to 1$. Our theory thus goes be-
yond previous, deterministic mean-field dynamics for first-
and second-order cumulants [22] or pseudocumulants [53]
that similarly captured the quenched variability of synaptic
weights but assumed the macroscopic limit $N \to \infty$.

The analysis of the new MF2 dynamics showed that the
quenched connectivity disorder for $p < 1$ has a drastic effect
on the neural variability, the response to time-dependent stim-
uli, and the stability of the network in biologically relevant
regimes. Specifically, our main findings with regard to the
effects of a finite connectivity $0 < p < 1$ and finite network
size $N < \infty$ are the following: (i) A finite variance in the input
currents ($\sigma^2 > 0$) caused by the quenched random connectiv-
ity with finite $C = pN$ leads to an effective broadening of the
neuron transfer function, but the stationary mean firing rate
remains rather unaffected [Figs. 3(a) and 8(b)]. (ii) In the pres-
ence of synaptic delays, the stability properties of fixed points
of the dynamics and locations of bifurcations can change
considerably. In particular, we showed that the connectivity
disorder can stabilize networks in regimes where an MF1
theory corresponding to a mean-connectivity approximation
would predict oscillatory dynamics [Figs. 3(b) and 3(c)]. (iii)
The population-rate response to high-frequency stimuli and
(iv) the variance of the population rate are significantly lower
than predicted by the MF1 theory [Figs. 3(d), 3(e) and 5].
Moreover, (v) in the presence of shared noise due to common
external input, an increase in the mean external input can
actually lead to a *decrease* of the variance of the population
rate (Figs. 1 and 6)—an effect which again is not captured in
the MF1 theory.

Overall, we provided simple analytical explanations for
the effects of finite connectivity and finite network size in
terms of the slope of the population transfer function and the
multiplicative character of the finite-size noise. In particular,
we found that a decreasing connection probability $p$ lowers
the slope $F_h$ through an increase of the membrane potential
variance $\sigma^2$ and that, at finite network size $N$, the neural vari-
ability has a rate-dependent contribution that is proportional
to the product of mean stimulus strength (or firing rate) and
the slope of the transfer function [Fig. 6(d) and Eq. (46)].

Of note, our mesoscopic approach differs from previous
theoretical studies that focused on how quenched random con-
nectivity shapes microscopic dynamics [48,72,75,76]. Neural

networks with quenched disorder can display a wide reper-
toire of dynamical states even under homogeneous random
connectivity [76], as considered here. In particular, inhibition-
dominated networks operating near the edge of chaos produce
rich single-neuron dynamics that may support complex com-
putations [76,77]. Importantly, this critical dynamical state
refers to the microscopic dynamics, as characterized by broad
distributions of pairwise correlations [76,77] and neuronal
timescales [78]. The annealed approximation underlying our
theory cannot reproduce these microscopic features. How-
ever, this limitation is not problematic for our mesoscopic
description of the coarse-grained population activity. As al-
ready noted in Ref. [76], the mesoscopic population activity is
largely unaffected by the critical dynamics at the microscopic
scale. We have furthermore shown here that, at the mesoscopic
scale, the annealed approximation captures surprisingly well
both the dynamics and their fluctuations across the full range
of connection probabilities.

To derive a simple mean-field model with fluctuations that
highlights the effects of finite connectivity and finite size,
we made several simplifications. In particular, we used (i) an
annealed approximation of the quenched random connectivity,
(ii) only one cell type (single population), and (iii) a first-order
Poisson model of the spiking-neuron dynamics. In the follow-
ing, we will briefly discuss these simplifications and possible
extensions.

In our "dynamically annealed" approximation, we assumed
that each presynaptic spike would be randomly distributed
among all neurons according to the connection probabil-
ity $p$, independent of the actual, fixed (quenched) synaptic
connectivity matrix of a given network. A systematic com-
parison of our theory with simulations of quenched networks
showed that our MF2 theory captures reliably the behavior
of quenched networks with respect to the main findings de-
scribed above. Deviations of the annealed network (from the
quenched network) and, correspondingly, the MF2 dynam-
ics, are noticeable for several statistics, though. However,
these deviations were considerably smaller than the respective
deviations of the mean-connectivity network and the corre-
sponding MF1 dynamics used in previous studies [36,47,79].
The remaining deviations of the annealed network can most
likely be attributed to the neglect of temporal correlations in
the recurrent input caused by the annealing, i.e., the incessant
resampling of the connections in time. This explanation would
be consistent with the underestimation of the low-frequency
power by the annealed network which we observed in Fig. 4.
For infinitely large, sparsely connected networks, a proper
account of temporal correlations can, in principle, be obtained
through a self-consistent treatment of the autocorrelations of
the recurrent fluctuations [62,75], which typically adds low-
frequency power [62]. Such self-consistent theory would lead
to colored noise rather than white noise for the incoherent
fluctuations in Eq. (12). However, how to build a stochastic,
dynamical population model for finite-size, random networks
that accounts for temporal correlations in a self-consistent
manner is an open theoretical problem which goes beyond the
scope of the present study.

Interestingly, perhaps not surprisingly, the comparison of
our theory with simulations of "dynamically annealed" net-
works shows that it can be considered almost exact for the
annealed case. One may ask whether a biological realization of such a “dynamically annealed” network, where our theory becomes exact, might actually exist. A possible origin of such randomly distributed spikes might be the probabilistic vesicle release: A fully connected network where synapses would transmit a spike with probability p would be an exact realization of the system we describe by our approximation. Thus, the annealed network model may also be regarded as a simple caricature of a spiking neural network with probabilistic synaptic transmission. Another type of annealing occurs in networks with full connectivity but dynamic synapses (short-term synaptic plasticity). Recently, such dynamically changing networks have also been reduced to low-dimensional Langevin dynamics based on a similar Poisson neuron model [58,59].

For the sake of a systematic and thorough development of the theory, we here began with the case of one population corresponding to a single neuron type; especially, we focused on the case of an inhibitory network. However, our theory suggests a straightforward extension to multiple populations. Cortical networks typically consist of several cell types and are generally described as excitatory-inhibitory networks that are, in the simplest case, organized as two coupled populations. Rate models of such two-population E-I networks remain a popular tool to describe network dynamics up to this day, whether to characterize inhibition stabilization in cortical networks [80] or spatial patterns of oscillatory activity observed in motor cortex [79,81]. It will certainly be worthwhile to investigate how the combined effects of random connectivity and network size we began to explore here may shape the more complicated dynamics of networks that comprise more than one population. Furthermore, a multipopulation extension would offer the possibility to model specific cortical microcircuits, such as canonical circuit models [82] or cortical column models with multiple interneuron types [41], on the mesoscopic scale, which accounts for fluctuations.

To develop the theory, we used Poisson neurons based on a one-dimensional dynamics for the membrane potential of each neuron. While a Poisson model neglects spike-history effects such as neuronal refractoriness, it allowed us to derive a low-dimensional, closed set of mean-field equations for the first- and second-order cumulants. We expect that in regimes where the neurons firing is Poissonian, similar observations can be made in more biophysically grounded models of spiking networks. While this is beyond the scope of this work, it would certainly be worthwhile to systematically investigate the effects described here in simulations of such networks.

The main difference of our Poisson model compared to more realistic integrate-and-fire (IF) models is the absence of a reset mechanism (in fact, with reset, our model would be equivalent to a leaky IF model with escape noise [28,31,33,83,84]). While the resetting in IF neurons captures refractoriness, it significantly complicates the population dynamics. Specifically, the population dynamics of IF neurons is usually infinite dimensional, represented by partial differential equations [29,84], integral equations [23], infinite systems of SDEs [66,85], or, in the case of finite N, stochastic versions of these [33,36,66,86]. In principle, it seems possible to carry over our theory to these population equations by again using the same annealed approximation. It would be interesting to

<figure><img src="image_10.png" /><figcaption>FIG. 7. Variance of $r$ and $\bar{h}$ for a network of LIF neurons with a rectified linear hazard function $\phi(h) = \beta_{\text{lin}}[h - \theta_{\text{lin}}]_+$ and a post-spike reset to $h_{\text{reset}}$ as a function of connection probability $p$. The LIF network behaves qualitatively similar to the Poisson network without reset in this paper (cf. Fig. 5). Variances are calculated from network simulations of quenched (green) and annealed (red) networks for (a) $N = 1000$ and (b) $C = 100$ for different connection probabilities $p$. The rectified linear hazard rate was matched both by value and slope of the single neuron transfer function Eq. (2) at a reference point $h_{\text{ref}} = -0.25$ mV. Parameters: $h_{\text{reset}} = -15$ mV, $\beta_{\text{lin}} \approx 91.32 \text{ mV}^{-1}$, $\theta_{\text{lin}} \approx -0.37$ mV, $w = -1$ mV s, $\mu_0 = 10$ mV.</figcaption></figure>

study how well the annealed approximation works for neurons with refractoriness. One may speculate that similar second-order mean-field theories may provide a quantitatively correct account for networks, e.g., of generalized integrate-and-fire neurons in the presence of finite connectivity where effective first-order mean-field theories have already proven to be quantitatively precise for fully connected networks [33,36,79]. While a systematic analysis of the effect of resets is beyond the scope of this paper, we present here some first simulation results. We consider a quenched network, in which the input potential is reset to $h_{\text{reset}}$ after each spike. With reset, the single-neuron hazard rate no longer requires an upper bound, and we therefore adopt a rectified linear hazard rate [28]. As an example, we repeat the simulations from Fig. 5 with reset. The results (Fig. 7) show that the variances of the population rate $r$ and mean input potential $\bar{h}$ as a function of the connection probability $p$ behave qualitatively similarly to the case without reset (Fig. 5). Specifically, for fixed $N$, the variance of $r$ increases with $p$, while the variance of $\bar{h}$ decreases. For fixed $C$, the variance of $r$ increases, and the variance of $\bar{h}$ exhibits a maximum. Importantly, the annealed approximation now matches the quenched network much more closely than in the absence of reset. This aligns with our earlier explanation that discrepancies between the quenched network and the annealed approximation arise from neglecting temporal correlations in the incoherent recurrent input due to annealing. The reset reduces the enhanced low-frequency power of the model without reset, effectively whitening the spike-train power spectra and thus reducing temporal correlations in the input. In other words, there are less temporal correlations to be removed by the annealing, leading to a better fit.
For simulations and concreteness, we made some further specific modeling choices: For example, we studied random networks with fixed in-degree, i.e., an identical number of presynaptic connections, and chose a specific transfer function (or f-I curve) $\phi(h)$ that relates the instantaneous firing rate to the intensity or input current $h$. While our specific choice of the transfer function also enabled us to analytically evaluate the effect of a finite variance in the input currents on the population rate, the reduction of the spiking network dynamics to a low-dimensional system of coupled SDEs should remain valid for arbitrary choices of the transfer function, albeit probably more difficult to compute. The assumption of fixed in-degree allowed us to assume that the mean input current as well as the strengths of the fluctuating parts are identical among neurons, but in practice we observed that a purely random Erdős-Rényi network shows almost identical statistics as those described in the previous sections for networks with fixed in-degree.

Here we showed how to derive low-dimensional stochastic models that describe fluctuations in a complex, biological model system, namely a finite-size spiking neural networks with random connectivity. Already in the simple case of one population and Poissonian spiking, we found highly nontrivial effects such as the nonmonotonic dependence of the variance of the recurrent input on connection probability [Fig. 5(b)] and the weak suppression of variability by external stimuli [Fig. 6(b)]. Thus, even though a single population and Poisson spiking dynamics may be too simple to quantitatively model certain phenomena in cortex (such as a stronger suppression of variability [61]), the fact that our MF2 theory explains nontrivial fluctuation effects demonstrates the power of the approach as a proof of principle. Importantly, with the biological extensions discussed above, we believe that our theory may pave the way for a useful modeling framework for variability in cortical circuits at the mesoscopic population level.

## ACKNOWLEDGMENT

We are grateful to Jakob Stubenrauch for useful comments on the manuscript.

## DATA AVAILABILITY

The data that support the findings of this article are openly available [87].

## APPENDIX A: MESOSCOPIC DYNAMICS IN THE PRESENCE OF HETEROGENEOUS EXTERNAL DRIVE

An additional source of disorder in the microscopic models is the heterogeneity in the parameters. For Eq. (1), we include heterogeneity in the external drive as follows:

$$ \mu_{\text{ext},i}(t) = \mu(t) + \hat{\mu}_i. \quad (A1) $$

The random variables $\hat{\mu}_i$ for each neuron $i$ have mean zero and are independently and identically distributed with probability density $\rho(\hat{\mu})$. They can be interpreted as heterogeneity of the external drives, the resting potentials or the thresholds $\vartheta$ of the different neurons. Importantly, the variables $\hat{\mu}_i$ are constant in time, and thus represent quenched disorder. The derivation of the microscopic annealed network (12) does not explicitly depend on the value of $\mu(t)$ and we therefore obtain the same results with $\mu(t)$ replaced with Eq. (A1). Specifically, Eq. (15)

is modified as follows:

$$ \tau \frac{dh_i}{dt} = -h_i + \hat{\mu}_i + f(t) + g(t)\zeta_i(t). $$

Splitting off the heterogeneous part from $h_i(t)$ yields a new variable $x_i(t) = h_i(t) - \hat{\mu}_i$ that obeys

$$ \tau \dot{x}_i = -x_i + f(t) + g(t) \zeta_i(t). $$

The last equation is of the same form as Eq. (15) and can be dealt with accordingly. For the self-consistent closure, the probability density of the input potentials $h_i$ is needed for the calculation of $\langle \phi(h_i(t)) \rangle$. At each time point we assume $h_i(t) = x_i(t) + \hat{\mu}_i$ is the sum of two independent random variables and thus its density is the convolution $g_{\bar{x}, \sigma_x^2} * \rho$. For simplicity, we assume that the heterogeneity is Gaussian distributed $\hat{\mu}_i \sim N(0, \sigma_\mu^2)$ with a given fixed variance $\sigma_\mu^2$. In this case, we have a Gaussian probability for the distribution $h_i$ with mean $\bar{x}(t)$ and variance $\sigma_x^2 + \sigma_\mu^2$. We therefore have

$$ \begin{aligned} \tau \frac{d\bar{x}}{dt} &= -\bar{x} + \mu(t) + w \left[ r(t-d) + \sqrt{\frac{r(t-d)}{N}} \eta(t) \right] \\ \tau \frac{d\sigma_x^2}{dt} &= -2\sigma_x^2 + \frac{w^2(1-p)}{\tau p N} r(t-d) \\ \tau \frac{d\xi}{dt} &= -\xi + \sqrt{2\tau G(\bar{x}(t), \sigma^2(t))} \zeta(t), \end{aligned} $$

where

$$ \begin{aligned} r(t) &= F(\bar{x}(t), \sigma^2(t)) + \frac{1}{\sqrt{N}} \xi(t) \\ \sigma^2(t) &= \sigma_x^2(t) + \sigma_\mu^2. \end{aligned} $$

In other words, we obtain the very same mean-field dynamics as in the case without heterogeneity, Eq. (20), but with an increased variance $\sigma^2$. See also the classical mean-field theory of Amari [22] for a similar treatment of heterogeneity.

## APPENDIX B: FIXED-POINT SOLUTIONS

At equilibrium, the sparse limit, Eq. (22), gives

$$ \bar{h}_0 = \mu_0 + w r_0, \quad \sigma_0^2 = \frac{w^2}{2\tau C} r_0, \quad r_0 = F(\bar{h}_0, \sigma_0^2), \quad (B1) $$

which yields the fixed-point equation (25) for $h_0$.

Because we assume an inhibitory network ($w < 0$), the linear function on the left-hand side of Eq. (25) has a negative slope. If we choose $w$ such that the intersection occurs for values $\bar{h}_0$ lower than the inflection point $\vartheta = 0$ of the hazard function $\phi(h)$ [corresponding to the biologically relevant convex part of $\phi(h)$], then from the graphical solution it is clear that the second-order MF dynamics exhibits a fixed-point with lower mean input potential $\bar{h}_0$ and higher population firing rate $r_0$ compared to the first-order MF dynamics. This effect becomes stronger for smaller connection probability $p$.

The mean stationary firing rates of the quenched and annealed microscopic models are well captured by the numerical solution of the fixed-point equation (25) [Figs. 8(a)–8(c)]. As a function of the external drive $\mu_0$, we observe a strikingly linear dependence between the saturation regimes at low and high values of $\mu_0$ [Fig. 8(a)]. Furthermore, we find that the firing rate is roughly constant when the connection probability
<figure><img src="image_11.png" /><figcaption>FIG. 8. Mean stationary population rate of the quenched and annealed networks and fixed-point solutions of the noiseless mean-field dynamics. Population averaged firing rates as a function of (a) the external input $\mu_0$, (b) the connection probability $p$ for $C = 100$ constant, and (c) the coupling strength $w$. Orange (blue) line: MF1 (MF2); brown line: MFsparselim; red squares (green circles); annealed (quenched) network simulation; dotted line: analytical approximation (26). The inset in (c) shows the same data, but for a larger range of $w$ in a double logarithmic scaling. Gray dashed line: second-order theory $r_\infty$ for the firing rate in the limit $w \to -\infty$ [Eq. (B5)]. Note: MFsparselim is identical to MF2 in (a) and (c). Parameters: $N = 1000$, $\beta = 5 \text{ mV}^{-1}$ (a) $p = 0.1$ and $w = -1 \text{ mV s}$; (b) $\mu_0 = 10 \text{ mV}$, $w = -1 \text{ mV s}$; and (c) $p = 0.1$, $\mu_0 = 10 \text{ mV}$.</figcaption></figure>

$p$ is varied while $N$ is fixed [Fig. 8(b)]. The weak dependence on $p$ is slightly better predicted by the second-order MF theory than by the first-order MF theory, which, by definition, has no dependence on $p$. We also observe a power-law behavior of the mean firing rate in both microscopic simulations and fixed-point solutions if the coupling strength is not too strong [Fig. 8(c)]. We note that, in general, the relative error of the first-order MF dynamics (and hence the mean-connectivity network) with respect to the stationary mean population rate is rather small compared to the much larger error with respect to the second-order statistics reported in Fig. 1.

## 1. Analytical approximation of fixed points

The linear dependence on $\mu_0$, the weak dependence on $p$ and the power-law dependence on $w$ can be explained with our MF theory as follows. First, we notice that for a sufficiently large steepness $\beta \gg 1/\sigma_0^2$ of the single-neuron transfer function $\phi$, the second-order MF dynamics remains practically unaltered if we take the limit $\beta \to \infty$. In this limit, the single-neuron transfer function tends to $r_m\theta(\bar{h})$, where $\theta$ denotes the Heaviside step function. Indeed, for the values of $\beta$ considered in this paper, taking the limit $\beta \to \infty$ does not noticeably alter the effective nonlinearity $F$ and hence the second-order MF dynamics. The limiting value of the fixed point can be calculated analytically from the graphical representation [Fig. 3(a)] which yields an approximate analytical expression for the firing rate in the first-order MF dynamics,

$$r_0 \approx \begin{cases} 0, & \mu_0 \le 0 \\ -\frac{\mu_0}{w}, & 0 < \mu_0 < -wr_m \\ r_m, & \mu_0 \ge -wr_m. \end{cases} \quad (B2)$$

This analytical expression predicts that $r_0$ is inversely proportional to the coupling strength in the balanced regime $0 < \mu_0 < -wr_m$. Indeed, the numerical fixed-point solution for the first-order MF dynamics exhibits an excellent quantitative agreement with this prediction [Fig. 8(c)].

Second, we note from the geometry of the graphical solution of the fixed points that the intersection point of the first-order MF dynamics serves as a rough approximation for the firing rate $r_0$ of the second-order MF dynamics as well

[Fig. 3(a)] unless $w$ is extremely negative as discussed below. Indeed, the horizontal stretching of the transfer function in the second-order theory mainly affects the value of the fixed point $\bar{h}_0$ rather than the firing rate $r_0$. More precisely, a sufficient condition for the validity of our approximation is that $-\mu_0/w \gg (1-p)/(\tau C)$ and $-\mu_0/w > 0.023r_m$ (see Appendix B 3). We expect that the approximation (26) slightly underestimates the firing rates if the intersection is below the inflection point ($\bar{h}_0 < 0$) and slightly overestimates the firing rates if the intersection is above the inflection point ($\bar{h}_0 > 0$). Our approximation Eq. (26) well explains the piecewise linear behavior of the firing rate as a function of the external input as observed in Fig. 8(a). Furthermore, Eq. (26) also explains why the firing rates are roughly independent of the connection probability $p$ [Fig. 8(b)], and why we see the power-law $w^{-1}$ even in nonfully connected networks [Fig. 8(c)] unless $w$ is extremely negative (this case is treated below). In general we note that, compared to the first-order MF dynamics, the second-order MF dynamics shows a better agreement with the firing rate in microscopic simulations; in particular, it fits the annealed network almost perfectly.

Because closed-form analytical approximation of the fixed points will be crucial for a qualitative—if not quantitative— discussion of the linear analysis around fixed points in the following sections, we provide here also approximations for the fixed-point values $\bar{h}_0$ and $\sigma_0^2$. Using our approximation for the firing rate, Eq. (26), we obtain from Eq. (B1) an approximation for the dispersion:

$$\sigma_0^2 \approx -\frac{w\mu_0(1-p)}{2\tau Np} \quad (B3)$$

if $0 < \mu_0 < -wr_m$. To obtain a corresponding approximation for $\bar{h}_0$ in the second-order MF theory, we solve the fixed-point equation $r_0 = F(\bar{h}_0, \sigma_0^2)$ [cf. Eq. (B1)] for $\bar{h}_0$ using our approximations for $r_0$ and $\sigma_0^2$ resulting in

$$\bar{h}_0 \approx \Phi^{-1}\left(-\frac{\mu_0}{wr_m}\right) \sqrt{\beta^{-2} + \sigma_0^2}. \quad (B4)$$

## 2. Limit of strongly inhibitory networks

Finally, the regime of strong inhibitory coupling reveals a qualitative difference between the first- and second-order
mean-field prediction. In the limit $w \to -\infty$, the right-hand side of Eq. (25) is asymptotic to

$$r_m \Phi \left( -\sqrt{\frac{2\tau pN(\bar{h}_0 - \mu_0)}{(1-p)w}} \right),$$

and thus the mean stationary firing rate in a strongly inhibitory network, is given by the solution $r_\infty$ of the limiting equation,

$$r_\infty = r_m \Phi \left( -\sqrt{\frac{2\tau pN r_\infty}{1-p}} \right). \quad (B5)$$

The graphical solution of this equation clearly shows a vanishing firing rate in the first-order MF theory ($p \to 1$), whereas the limiting equation has a nonvanishing solution for $p < 1$ [second-order MF theory, Fig. 8(c), inset]. The nonvanishing, limiting firing rate in the second-order MF theory is in good agreement with simulations of the annealed and quenched microscopic model. In contrast, the first-order MF theory generally underestimates the firing rates in strongly inhibitory networks: For $w \to -\infty$, the first-order MF dynamics has a rate that converges to zero like $-\mu_0/w$ as predicted by Eq. (26).

### 3. Condition for the validity of the firing-rate approximation (26)

In Eq. (26), a simple approximation for the fixed-point firing rate is given, namely $r_0 \approx -\mu_0/w$, which we used for an inhibitory network ($w < 0$) in the nonsaturated regime ($0 < \mu_0 \le -wr_m$). This approximation corresponds to the solution of the first-order MF theory in the limit $\beta \to \infty$. We argued that under certain sufficient conditions, the approximation is good for the second-order MF dynamics as well. Here we provide a justification for these conditions. For biological reasons, we make a slightly stronger assumption on the working point of our system: We assume that the fixed point is in the convex part of the transfer function, i.e., $h_0 < 0$ and $r_0 < r_m/2$, as drawn in Fig. 3(a).

For the firing-rate approximation to be applicable also to the second-order MF dynamics, we require that the relative difference between the firing rates of the first- and second-order MF theory is small:

$$r_0^{2\text{nd}} - r_0^{1\text{st}} \ll r_0^{1\text{st}}. \quad (B6)$$

Here $r_0^{1\text{st}} = -\mu_0/w$ and $r_0^{2\text{nd}} = F(h_0, \sigma_0^2)$. Furthermore, we have used that $r_0^{2\text{nd}} > r_0^{1\text{st}}$ because of our assumption, $h_0 < 0$ [Fig. 3(a)]. To estimate the left-hand side of (B6) from above, we want to obtain an upper bound approximation of the rate $r_0^{2\text{nd}} = (h_0 - \mu_0)/w$. To this end, we lower bound the input potential $h_0$ (note that $w < 0$). Because $h_0$ is expected to be on the order of the width $\sigma_0$ of the transfer function, we write the lower bound as $-a\sigma_0 < h_0$, where $a$ is a positive number that is still unknown. This lower bound of $h_0$ yields an upper bound estimate for the mean firing rate: $r_0^{2\text{nd}} < (-a\sigma_0 - \mu_0)/w$. Hence, the condition (B6) is surely met if

$$\frac{-a\sigma_0 - \mu_0}{w} - \frac{-\mu_0}{w} = -a \frac{\sigma_0}{w} \ll -\frac{\mu_0}{w},$$

hence $a\sigma_0 \ll \mu_0$. Unfortunately, we do not know the exact value of $\sigma_0$ for the second-order model. However, if our approximation holds true, then we can self-consistently use our approximation for $\sigma_0^2$, Eq. (28), and the condition becomes

$$a \sqrt{\frac{-w(1-p)}{2\tau \mu_0 C}} \ll 1. \quad (B7)$$

The larger the value of $a$, the more conservative this sufficient condition becomes. However, to obtain a less conservative condition, we could lower the value of $a$, as long as $h_0 > -a\sigma_0$. We check for this latter condition in the same manner as above by self-consistently using the approximation of $h_0$ in Eq. (27):

$$h_0 \approx \Phi^{-1}\left(\frac{-\mu_0}{w r_m}\right) \sigma_0 > -a\sigma_0,$$

where we have used the limit $\beta \to \infty$. Eventually, this leads to

$$\frac{-\mu_0}{w r_m} > \Phi(-a). \quad (B8)$$

The last step was possible, because $\Phi$ is strictly monotonously increasing. If there exists $a > 0$ such that both Eqs. (B7) and (B8) are fulfilled, then we consider the approximations for the fixed points to be valid. For the parameter choice used in our simulations, $a = 2$ was sufficient. The numerical evaluation of the error function yields $\Phi(-2) \approx 0.023$. Fixing $a = 2$ yields the sufficient condition reported in Sec. IV A.

## APPENDIX C: ESTIMATION OF THE RATE VARIANCE FROM POPULATION ACTIVITIES

The variance of the population firing rate is not directly accessible from measurements in real biological networks. However, we can estimate the rate variance by the variance of the empirical population activity. We assume that in small time bins $(t, t + \Delta t)$, neurons fire spikes independently with conditional intensities $\lambda_i(t)$ given the past. Then, from the law of total variance and the conditionally Poisson statistics, we have for the variance of the total number of spikes in a time bin

$$\begin{aligned} \text{var}(\Delta Z) &= \langle \text{var}(\Delta Z | \{\lambda_i\}) \rangle + \text{var}(\langle \Delta Z | \{\lambda_i\} \rangle) \\ &= N \langle r \rangle \Delta t + N^2 \text{var}(r) \Delta t^2, \end{aligned}$$

where $r(t) = N^{-1} \sum_{i=1}^N \lambda_i(t)$ is the population firing rate and $\langle \cdot \rangle$ denotes the trial average. Hence,

$$\begin{aligned} \text{var}(r) &= \frac{\text{var}(\Delta Z)}{N^2 \Delta t^2} - \frac{\langle r \rangle}{N \Delta t} \\ &= \text{var}(A_N(t, \Delta t)) - \frac{\langle A_N(t, \Delta t) \rangle}{N \Delta t}, \end{aligned}$$

where we used the definition of the empirical population activity, Eq. (3). For stationary data, the trial average can be replaced by the time average.
## APPENDIX D: PARTIAL DERIVATIVES OF THE FUNCTION F AND APPROXIMATIONS THEREOF

Using the definition of the function F defined in Eq. (23), we find for the partial derivative with respect to h at the fixed point $(\bar{h}_0, \sigma_0^2)$,

$$F_h := \frac{\partial F}{\partial h}(\bar{h}_0, \sigma_0^2 + \sigma_\mu^2) = \frac{r_m \beta \exp\left[-\frac{\beta^2 h_0^2}{2(1+\beta^2(\sigma_0^2+\sigma_\mu^2))}\right]}{\sqrt{2\pi(1+\beta^2(\sigma_0^2+\sigma_\mu^2))}}$$

Note that $\sigma_\mu^2 = 0$ in the absence of heterogeneity of $\mu_i$. Similarly, for the partial derivative with respect to the second argument, we find

$$F_\sigma := \frac{\partial F}{\partial \sigma^2}(\bar{h}_0, \sigma_0^2 + \sigma_\mu^2) = \frac{-r_m \beta^3 h_0 \exp\left[-\frac{\beta^2 h_0^2}{2(1+\beta^2(\sigma_0^2+\sigma_\mu^2))}\right]}{2\sqrt{2\pi}\sqrt{1+\beta^2(\sigma_0^2+\sigma_\mu^2)^{-3}}}$$

The approximate expressions Eqs. (27) and (28) for $h_0$ and $\sigma_0^2$ at the fixed point can be used for further simplifications of the partial derivatives:

$$F_h \approx \frac{r_m \beta Q(-\frac{\mu_0}{r_m w})}{\sqrt{2\pi} \sqrt{1 + \beta^2(\sigma_0^2 + \sigma_\mu^2)}}$$

where $\sigma_0^2$ and the function $Q$ are given by Eqs. (28) and (38), respectively. Analogously, we can approximate the partial derivative with respect to $\sigma^2$ at the fixed point:

$$F_\sigma \approx \frac{-r_m \beta^2 \Phi^{-1} \left(\frac{-\mu_0}{r_m w}\right) Q \left(\frac{-\mu_0}{r_m w}\right)}{2\sqrt{2\pi} (1 + \beta^2(\sigma_0^2 + \sigma_\mu^2))}$$

## APPENDIX E: LINEAR STABILITY ANALYSIS

For the linear stability analysis, we consider the linearized system with no noise and constant $\mu(t) = \mu$ which reads

$$\tau \frac{d}{dt} \delta h(t) = -\delta h(t) + w \delta r(t - d), \quad (E1a)$$

$$\tau \frac{d}{dt} \delta \sigma^2(t) = -2 \delta \sigma^2(t) + \frac{w^2(1-p)}{\tau C} \delta r(t - d), \quad (E1b)$$

$$\delta r(t) = F_h \delta h(t) + F_\sigma \delta \sigma^2(t), \quad (E1c)$$

where we neglected the variable $\xi$ as it decays exponentially. To solve the linearized system, we insert the exponential ansatz

$$\delta h(t) = \hat{h}(\lambda)e^{\lambda t} \\ \delta \sigma^2(t) = \hat{\sigma}(\lambda)e^{\lambda t},$$

with a constant parameter $\lambda$. Solving Eq. (E1b) for $\hat{\sigma}$ yields

$$\hat{\sigma} = \frac{w^2 B F_h e^{-\lambda d}}{\lambda \tau + 2 - w^2 B F_\sigma e^{-\lambda d}} \hat{h}, \quad B := \frac{1-p}{\tau C}. \quad (E2)$$

We insert this expression for $\hat{\sigma}$ into the first equation (E1a) and obtain an equation for the eigenvalue

$$\lambda \tau = -1 + w F_h e^{-\lambda d} + \frac{\hat{F}_\sigma w F_h e^{-2\lambda d}}{2 + \lambda \tau - \hat{F}_\sigma e^{-\lambda d}}, \quad (E3)$$

where we used the abbreviation $\hat{F}_\sigma := w^2 B F_\sigma$. This equation depends on the fixed points through the functions $F_h$ and $F_\sigma$. To find an oscillatory instability, we look for a Hopf bifurcation, at which $\lambda = i\omega$ for some real-valued frequency $\omega$. Applying the condition on the complex-valued equation (E3) we obtain two conditions, one for the amplitude,

$$1 = \frac{4 + \omega^2 \tau^2}{1 + \omega^2 \tau^2} \frac{w^2 F_h^2}{(2 - \cos(\omega d) \hat{F}_\sigma)^2 + (\omega \tau + \sin(\omega d) \hat{F}_\sigma)^2}, \quad (E4)$$

and one for the phase,

$$\arctan(\omega \tau) = \arctan\left(\frac{\omega \tau}{2}\right) + \arg(w F_h) - \omega d + 2\pi k - \operatorname{atan2}(2 - \hat{F}_\sigma \cos(\omega d), \omega \tau + \hat{F}_\sigma \sin(\omega d)).$$

Both must be satisfied simultaneously. Here atan2 is the two-argument arctangent and $k \in \mathbb{N}_0$. In the $(d, w)$-parameter space, these conditions correspond to a family of curves, one curve for each $k$. At the $k$th curve an oscillatory perturbation with frequency $\omega_k$ becomes unstable. For the overall instability boundary it is sufficient that at least one mode becomes unstable. We found empirically that the overall instability boundary is given by the mode $k=0$, whereas higher modes yield regions of instability that are contained already in the instability region given by $k=0$. For the solution of both conditions a numerical solution of the fixed-point equations is necessary. In Fig. 3(b) we scan through each parameter point and solve the phase condition for $\omega$. In the investigated parameter range we were always able to find a unique solution for the phase condition for each mode. We then determined the boundary in the $(d, w)$-space that separated the regions where the right-hand side of the amplitude condition, Eq. (E4), is larger and smaller than 1, respectively. The equations are simpler for the first-order MF dynamics where $\hat{F}_\sigma = 0$. The simplified amplitude condition reads in this case

$$1 = \frac{w^2 F_h^2}{1 + \omega^2 \tau^2} \Leftrightarrow \omega = \pm \frac{\sqrt{w^2 F_h^2 - 1}}{\tau}.$$

This explicit solution for $\omega$ can be used in the phase condition to determine the value of the delay at instability boundary for a given coupling strength $w$:

$$d = \frac{1}{\omega}(-\arctan(\omega\tau) + \arg(wF_h) + 2\pi k).$$

## APPENDIX F: SUSCEPTIBILITY MATRIX

Here we present explicit expressions for the susceptibility matrix $\tilde{\chi}$ defined in Eq. (33). Using Cramer's rule we obtain

$$\tilde{\chi}_{11} = \frac{1}{D} \left( i\omega + \frac{2}{\tau} - F_\sigma \tilde{\beta} e^{-i\omega d} \right), \quad (F1a)$$

$$\tilde{\chi}_{12} = \frac{1}{D} F_\sigma \tilde{\alpha} e^{-i\omega d}, \quad (F1b)$$

$$\tilde{\chi}_{13} = \frac{1}{D} \tilde{\gamma} \tilde{\alpha} e^{-i\omega d} \frac{i\omega + \frac{2}{\tau}}{i\omega + \frac{1}{\tau}}, \quad (F1c)$$

$$\tilde{\chi}_{21} = \frac{1}{D} F_h \tilde{\beta} e^{-i\omega d}, \quad (F1d)$$
$$ \tilde{\chi}_{22} = \frac{1}{D} \left( i\omega + \frac{1}{\tau} - F_h \tilde{\alpha} e^{-i\omega d} \right), \quad (\text{F1e}) $$

$$ \tilde{\chi}_{23} = \frac{1}{D} \tilde{\gamma} \tilde{\beta} e^{-i\omega d}, \quad (\text{F1f}) $$

$$ \tilde{\chi}_{31} = \tilde{\chi}_{32} = 0, \quad (\text{F1g}) $$

$$ \tilde{\chi}_{33} = \frac{1}{D} \left( i\omega + \frac{2}{\tau} \right) - \frac{1}{D} F_{\sigma} \tilde{\beta} e^{-i\omega d}, \quad (\text{F1h}) $$

$$ -\frac{i\omega + \frac{2}{\tau}}{D i\omega + \frac{1}{\tau}} F_h \tilde{\alpha} e^{-i\omega d}, \quad (\text{F1i}) $$

with the abbreviations

$$ \begin{aligned} D(\omega) = & \left( i\omega + \frac{1}{\tau} \right) \left( i\omega + \frac{2}{\tau} \right) - \left( i\omega + \frac{2}{\tau} \right) F_h \tilde{\alpha} e^{-i\omega d} \\ & - \left( i\omega + \frac{1}{\tau} \right) F_{\sigma} \tilde{\beta} e^{-i\omega d} \\ \tilde{\alpha} = & \frac{w}{\tau}, \quad \tilde{\beta} = \frac{w^2}{\tau^2} \left( \frac{1}{C} - \frac{1}{N} \right), \quad \tilde{\gamma} = \frac{1}{\sqrt{N}}. \end{aligned} $$

## APPENDIX G: VARIANCE AND POWER SPECTRAL DENSITY IN THE STATIONARY STATE

With the susceptibility we can express the power spectral density of the system. From Eq. (32), we obtain

$$ \langle \tilde{X}(\omega) \tilde{X}^*(\omega) \rangle = \tilde{\chi} \left\langle \left(B \tilde{\zeta} + \tilde{M}\right) \left(\tilde{\zeta}^* B^* + \tilde{M}^*\right) \right\rangle \tilde{\chi}^* $$

and with the definition of the power spectral density matrix $S(\omega)$, Eq. (39),

$$ S(\omega) = \tilde{\chi}(\omega)[S_M(\omega) + BB^*]\tilde{\chi}^*(\omega). $$

Here $S_M$ is the 3×3 spectral density matrix of the external stimulus. In our case, the external stimulus only acts on the $\tilde{h}$ variable and thus ($S_M$)$_{ij} = S_{\mu\mu}\delta_{1i}\delta_{1j}$ has only one nonzero entry with the power spectral density $S_{\mu\mu}(\omega)$ of the Gaussian stimulus $\mu(t)$. Explicitly,

$$ S_{ij}(\omega) = \frac{1}{\tau^2} \left[ S_{\mu\mu}(\omega) + \frac{w^2 r_0}{N} \right] \tilde{\chi}_{i1}(\omega) \tilde{\chi}_{j1}^*(\omega) + \frac{2G(\tilde{h}_0, \sigma_{h,0}^2)}{\tau N} \tilde{\chi}_{i3}(\omega) \tilde{\chi}_{j3}^*(\omega). $$

In the main text, we used this expression to compute the power spectrum of the population activity $A(t)$, given by Eq. (42). Taking the limit $p \to 1$ yields an explicit expression for the corresponding power spectrum of the first-order MF theory,

$$ S_{AA} \approx \frac{r_0}{N} + \frac{S_{\mu\mu} + w^2 \frac{r_0}{N} \phi_h^2 + 2w \frac{r_0}{N} \phi_h (1 - \cos(\omega d) w \phi_h)}{(\omega \tau - \sin(\omega d) w \phi_h)^2 + (1 - \cos(\omega d) w \phi_h)^2}. \quad (\text{G1}) $$

Here $\phi_h$ is the partial derivative $\partial\phi/\partial h$ at the fixed point.

In principle, the stationary variance of $h(t)$ can be computed from the integral over the power spectral density $S_{11}(\omega)$. However, the integral is difficult to solve for nonzero transmission delay and general spectral statistics $S_M(\omega)$ of the stimulus $\mu(t)$. If we restrict ourselves to the case of zero delay, $d=0$, and a Gaussian white noise stimulus,

$\mu(t) = \mu_0 + \sqrt{\tau \sigma_{\text{ext}}^2} \hat{\zeta}(t)$, then the linearized system can be rewritten as

$$ \frac{d}{dt} \mathbf{X} = \Gamma \mathbf{X} + \begin{pmatrix} \sqrt{\frac{\sigma_{\text{ext}}^2}{\tau}} \hat{\zeta}(t) \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} \frac{w}{\tau} \sqrt{\frac{r_0}{N}} \eta(t) \\ 0 \\ \sqrt{\frac{2G_0}{\tau}} \zeta(t) \end{pmatrix}, \quad (\text{G2}) $$

where $\Gamma = \text{T} + \text{W}$ or explicitly

$$ \Gamma = \begin{pmatrix} \frac{w}{\tau} F_h - \frac{1}{\tau} & \frac{w}{\tau} F_{\sigma} & \frac{w}{\tau} \frac{1}{\sqrt{N}} \\ \frac{w^2(1-p)}{\tau^2 N p} F_h & \frac{w^2(1-p)}{\tau^2 N p} F_{\sigma} - \frac{2}{\tau} & \frac{w^2(1-p)}{\tau^2 N p} \frac{1}{\sqrt{N}} \\ 0 & 0 & -\frac{1}{\tau} \end{pmatrix}. $$

The two Gaussian white noise processes in the first component of Eq. (G2) can be lumped together into a single Gaussian white noise process:

$$ \frac{d}{dt} \mathbf{X} = \Gamma \mathbf{X} + \hat{\mathbf{B}} \left( \eta(t) \right), \quad \hat{\mathbf{B}} = \begin{pmatrix} \sqrt{\frac{w^2 r_0}{\tau^2 N} + \frac{\sigma_{\text{ext}}^2}{\tau}} & 0 \\ 0 & 0 \\ 0 & \sqrt{\frac{2G_0}{\tau}} \end{pmatrix}. $$

This is a three-dimensional Ornstein-Uhlenbeck process for which the stationary variance is known to be the solution of the linear system of equations [88] (Lyapunov equation)

$$ \Gamma \sigma + \sigma \Gamma^T = -\hat{\mathbf{B}} \hat{\mathbf{B}}^T. \quad (\text{G3}) $$

Because the covariance matrix $\sigma$ is symmetric, $\sigma_{ij} = \sigma_{ji}$, we have a system of six linear equations for the variables $\sigma_{11}, \sigma_{12}, \sigma_{13}, \sigma_{22}, \sigma_{23}$, and $\sigma_{33}$. For a systematic solution, we notice that the matrix equation (G3) reads in the position (3, 3)

$$ 2(\Gamma_{31}\sigma_{13} + \Gamma_{32}\sigma_{23} + \Gamma_{33}\sigma_{33}) = -\frac{2G_0}{\tau}, $$

which reduces to $\sigma_{33} = G_0$ because $\Gamma_{31} = \Gamma_{32} = 0$ and $\Gamma_{33} = -1/\tau$. Therefore, we only need to solve a system of five linear equations

$$ \begin{pmatrix} 2\Gamma_{11} & 2\Gamma_{12} & 0 & 2\Gamma_{13} & 0 \\ 0 & 2\Gamma_{21} & 2\Gamma_{22} & 0 & 2\Gamma_{23} \\ \Gamma_{21} & \Gamma_{11} + \Gamma_{22} & \Gamma_{12} & \Gamma_{23} & \Gamma_{13} \\ 0 & 0 & 0 & \Gamma_{11} + \Gamma_{33} & \Gamma_{12} \\ 0 & 0 & 0 & \Gamma_{21} & \Gamma_{22} + \Gamma_{33} \end{pmatrix} \begin{pmatrix} \sigma_{11} \\ \sigma_{12} \\ \sigma_{22} \\ \sigma_{13} \\ \sigma_{23} \end{pmatrix} = \begin{pmatrix} \frac{w^2 r_0}{\tau^2 N} + \frac{\sigma_{\text{ext}}^2}{\tau} \\ 0 \\ 0 \\ -\frac{\Gamma_{13} G_0}{\tau} \\ -\Gamma_{23} G_0 \end{pmatrix}^T. $$

We can solve the last two equations separately as they only contain two variables:

$$ \left( \begin{array}{c} \sigma_{13} \\ \sigma_{23} \end{array} \right) = - \frac{G_0}{D} \left( \begin{array}{c} \Gamma_{13}(\Gamma_{22} + \Gamma_{33}) - \Gamma_{23}\Gamma_{12} \\ -\Gamma_{13}\Gamma_{21} + \Gamma_{23}(\Gamma_{11} + \Gamma_{33}) \end{array} \right), \quad (\text{G4}) $$

where $D = (\Gamma_{11} + \Gamma_{33})(\Gamma_{22} + \Gamma_{33}) - \Gamma_{21}\Gamma_{12}$. This equation gives us the values of $\sigma_{13}$ and $\sigma_{23}$. The remaining three equations of the form

$$ \begin{pmatrix} 2\Gamma_{11} & 2\Gamma_{12} & 0 \\ 0 & 2\Gamma_{21} & 2\Gamma_{22} \\ \Gamma_{21} & \Gamma_{11} + \Gamma_{22} & \Gamma_{12} \end{pmatrix} \begin{pmatrix} \sigma_{11} \\ \sigma_{12} \\ \sigma_{22} \end{pmatrix} = \begin{pmatrix} -\frac{w^2 r_0}{\tau^2 N} - \frac{\sigma_{\text{ext}}^2}{\tau} \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 2\Gamma_{13}\sigma_{13} \\ 2\Gamma_{23}\sigma_{23} \\ \Gamma_{23}\sigma_{13} + \Gamma_{13}\sigma_{23} \end{pmatrix} $$
have the solution

$$ (\sigma_{11}, \sigma_{12}, \sigma_{22})^T = \frac{-2}{\tau E} \left( \frac{w^2 r_0}{\tau N} + \sigma_{\text{ext}}^2 \right) \begin{pmatrix} \Gamma_{12}\Gamma_{21} - \Gamma_{22}(\Gamma_{11} + \Gamma_{22}) \\ \Gamma_{22}\Gamma_{21} - \Gamma_{21}^2 \end{pmatrix} \\ - \frac{4\sigma_{13}}{E} \begin{pmatrix} \Gamma_{13}[\Gamma_{21}\Gamma_{12} - \Gamma_{22}(\Gamma_{11} + \Gamma_{22})] + \Gamma_{23}\Gamma_{12}\Gamma_{22} \\ \Gamma_{13}\Gamma_{22}\Gamma_{21} - \Gamma_{23}\Gamma_{11}\Gamma_{22} - \Gamma_{13}\Gamma_{21}^2 + \Gamma_{23}\Gamma_{11}\Gamma_{21} \end{pmatrix} \\ - \frac{4\sigma_{23}}{E} \begin{pmatrix} -\Gamma_{23}\Gamma_{12}^2 + \Gamma_{13}\Gamma_{12}\Gamma_{22} \\ \Gamma_{23}\Gamma_{11}\Gamma_{12} - \Gamma_{13}\Gamma_{11}\Gamma_{22} - \Gamma_{23}[\Gamma_{11}(\Gamma_{11} + \Gamma_{22}) - \Gamma_{21}\Gamma_{12}] + \Gamma_{13}\Gamma_{11}\Gamma_{21} \end{pmatrix} \tag{G5} $$

with the determinant of the remaining 3×3 subsystem

$$ E = -4(\Gamma_{11}\Gamma_{22} - \Gamma_{12}\Gamma_{21})(\Gamma_{11} + \Gamma_{22}). $$

The formula for the variances of the mean input potential $\bar{h}$ and the firing rate $r$ can be more easily interpreted if some further, heuristic simplifications are made. For the parameters used in our study, we observed that $\Gamma_{11}$ is typically much larger than $\Gamma_{22}$, so that we could safely replace the terms $\Gamma_{11} + \Gamma_{22}$ in Eq. (G5) by $\Gamma_{11}$. Likewise, we found that the terms proportional to $\sigma_{13}$ and $\sigma_{23}$ in Eq. (G5) are small enough to be ignored. These simplifications yield

$$ \sigma_{11} \approx \frac{1}{2} \frac{1}{(1 - w F_h)} \left( \frac{w^2 r_0}{\tau N} + \sigma_{\text{ext}}^2 \right). $$

To simplify even further, we note that in our simulations we have that $|w|F_h \gg 1$, and hence

$$ \sigma_{11} \approx \left( \frac{-w r_0}{2 \tau N} + \frac{\sigma_{\text{ext}}^2}{-2w} \right) \frac{1}{F_h}. \tag{G6} $$

As mentioned in the main text, the variance of the population firing rate is dominated by the first term in Eq. (43), so that we obtain from Eq. (G6)

$$ \operatorname{var}(r) \approx \left( \frac{-w r_0}{2 \tau N} + \frac{\sigma_{\text{ext}}^2}{-2w} \right) F_h. \tag{G7} $$

The variance in the first-order MF theory can be either calculated directly from Eq. (20) or can be understood as a special case of the formulas Eqs. (G6) and (G7) for $p \to 1$ while keeping $N$ constant. The variances for the mean input and the population firing rate read

$$ \operatorname{var}(\bar{h}) = \frac{1}{2(1 - w\phi_h)} \left( \frac{w^2 r_0}{\tau N} + \sigma_{\text{ext}}^2 \right), \\ \operatorname{var}(r) = \frac{1}{2(1 - w\phi_h)} \left( \frac{w^2 r_0}{\tau N} + \sigma_{\text{ext}}^2 \right) \phi_h^2, $$

where $\phi_h$ is the derivative of the single neuron transfer function at the fixed point. For the approximation $1 - w\phi_h \approx -w\phi_h$ we arrive at the same form as Eqs. (G6) and (G7).

## APPENDIX H: CALCULATION OF THE FUNCTIONS F AND G

For the evaluation of Eqs. (23) and (24) for the given nonlinearity (2), we use two known formulas [89] for

Gaussian integrals,

$$ \int_{-\infty}^{\infty} \Phi(a + bx) g_{0,1}(x) dx = \Phi\left(\frac{a}{\sqrt{1 + b^2}}\right) \\ \int_{-\infty}^{\infty} \Phi^2(a + bx) g_{0,1}(x) dx = \Phi\left(\frac{a}{\sqrt{1 + b^2}}\right) \\ - 2T\left(\frac{a}{\sqrt{1 + b^2}}, \frac{1}{\sqrt{1 + 2b^2}}\right) $$

with $a, b \in \mathbb{R}$ and $g_{0,1}$ being the standard normal distribution. Here

$$ T(h, a) = g_{0,1}(h) \int_{0}^{a} \frac{g_{0,1}(hx)}{1 + x^2} dx $$

is the Owen's $T$ function. For our integral (23) we need to substitute $x := (h - \bar{h})/\sigma$

$$ F(\bar{h}, \sigma^2) = r_m \int_{-\infty}^{\infty} \frac{\Phi(\beta(h - \vartheta))}{\sqrt{2\pi}\sigma^2} \exp\left[-\frac{(h - \bar{h})^2}{2\sigma^2}\right] dh \\ = \frac{r_m}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \Phi(\beta(h - \vartheta) + \beta\sqrt{\sigma^2}x)e^{-x^2/2} dx \\ = r_m \Phi\left(\frac{\beta(h - \vartheta)}{\sqrt{1 + \beta^2\sigma^2}}\right). $$

The second integral (24) is computed analogously using the same substitution.

## APPENDIX I: THERE EXISTS EXACTLY ONE FIXED POINT FOR INHIBITORY NETWORKS WITH POSITIVE EXTERNAL DRIVE

For the calculation of the fixed points of the second-order model (20), we set the left-hand side to zero and solve the following equations for $h_0$:

$$ r_0 = F(h_0, \sigma^2(h_0)) = r_m \Phi\left(\frac{\beta h_0}{\sqrt{1 + \alpha(\mu_0 - h_0)}}\right) =: F(h_0), \tag{I1} $$

$$ r_0 = \frac{h_0 - \mu_0}{w}, \tag{I2} $$

where $\alpha = -\beta^2 w(1-p)/(2\tau C) > 0$. The strictly positive nonlinear function $F$ in Eq. (I1) can only intersect the linear function in Eq. (I2) at a value for the fixed point $h_0 < \mu_0$. Equation (I2) represents a strictly monotonically decreasing function of $h_0$ for $w < 0$. We will show that for $\mu_0 > 0$, the function $F(h)$ is strictly monotonically increasing for $h_0 < \mu_0$. The consequence is that both graphs will only intersect at
<table><caption>TABLE II. Model parameters and their values.</caption><thead><tr><th>Parameter</th><th>Value/relations</th><th>Description</th></tr></thead><tbody><tr><td>τ</td><td>20 ms</td><td>Time constant of the low-pass filter dynamics</td></tr><tr><td>r<sub>m</sub></td><td>100 Hz</td><td>Maximal firing rate</td></tr><tr><td>β</td><td>5 mV<sup>-1</sup></td><td>Steepness of the single-neuron transfer function (threshold sharpness)</td></tr><tr><td>ϑ</td><td>0 mV</td><td>Inflection point of single-neuron transfer function (soft threshold)</td></tr><tr><td>d</td><td>0 ms (unless otherwise noted)</td><td>Transmission delay</td></tr><tr><td>p</td><td>see figure captions</td><td>Connection probability</td></tr><tr><td>N</td><td>see figure captions</td><td>Number of neurons</td></tr><tr><td>C</td><td>C = pN</td><td>Number of incoming synapses</td></tr><tr><td>w</td><td>w = -1 mV s (unless noted otherwise)</td><td>Total coupling strength</td></tr><tr><td>J</td><td>J = w/C, (mVs)</td><td>Efficacy of single synapse</td></tr><tr><td>a<sub>ij</sub></td><td>(a<sub>ij</sub>) ~ Uniform({A ∈ {0, 1}<sup>N×N</sup>: Σ<sub>j=1</sub><sup>N</sup> A<sub>ij</sub> = C})</td><td>Adjacency matrix, a<sub>ij</sub> = 1[neuron j connected to i]</td></tr><tr><td>J<sub>ij</sub></td><td>J<sub>ij</sub> = J a<sub>ij</sub></td><td>Synaptic efficacy from neuron j to i</td></tr><tr><td>w<sub>ij</sub></td><td>w<sub>ij</sub> = NJ<sub>ij</sub>, (w<sub>ij</sub>) = w</td><td>Rescaled synaptic efficacy</td></tr></tbody></table>

exactly one point providing a unique fixed point. The proof for the monotonicity of $F(h_0)$ is straightforward. The derivative is

given by the derivative of the error function:

$$ \frac{dF}{dh} = \frac{r_m \beta \exp \left[ - \frac{\beta^2 h^2}{1 + \alpha(\mu_0 - h)} \right]}{\sqrt{2\pi} \sqrt{1 + \alpha(\mu_0 - h)}} \left( 1 + \frac{1}{2} \cdot \frac{\alpha h}{1 + \alpha(\mu_0 - h)} \right). $$

<table><caption>TABLE III. Performance summary of the mean-field theories for different parameter regimes. *Exact, apart from the minor diffusion approximation of the shot noise, which is valid for sufficiently large C and which is actually not really necessary: one could go back and replace the term $A(t) \approx r(t) + \sqrt{r(t)/N}\eta(t)$ by a Poissonian shot noise $A(t) = Pois(Nr(t)dt)/(Ndt)$, which would recover the exact dynamics.</caption><thead><tr><th>Parameter regime</th><th>MF1, Eq. (21)</th><th>MF2, Eq. (20)</th><th>MFsparselim, Eq. (22)</th></tr></thead><tbody><tr><td>N, C → ∞,<br>p = C/N = O(1)<br>(dense)</td><td colspan="3">τ (d/dt) h̄ = -h̄ + μ<sub>ext</sub>(t) + w r(t - d), r(t) = φ(h̄(t))<br>• Exact</td></tr><tr><td>N → ∞, C = O(1)<br>(sparse)</td><td>→ Eq. (J1)<br>✓ Stationary population rate<br>(for not too strong coupling)<br>✗ Time-dependent rate<br>response poor</td><td colspan="2">τ (d/dt) h̄ = -h̄ + μ<sub>ext</sub>(t) + w r(t - d)<br>τ (d/dt) σ<sup>2</sup> = -2σ<sup>2</sup> + w<sup>2</sup> / τC r(t - d),<br>r(t) = F(h̄(t), σ<sup>2</sup>(t))<br>• Exact for annealed network<br>• For quenched network:<br>✓ Stationary population rate better than MF1<br>✓ Synamical rate response</td></tr><tr><td>N < ∞, p = 1</td><td colspan="3">τ (d/dt) h̄ = -h̄ + μ<sub>ext</sub>(t) + w [r(t - d) + √(r(t - d)/N) η(t)]<br>r(t) = φ(h̄(t))<br>• Exact*</td></tr><tr><td>N < ∞, 0 < p < 1</td><td>→ Eq. (J2) with C = N<br>✓ No finite-size fluctuations<br>✗ Stationary and<br>time-dependent rate poor</td><td colspan="2">✓ Captures finite-size fluctuations<br>very well for annealed network<br>✓ Strong improvement for<br>quenched network compared<br>to MF1</td></tr></tbody></table>
The first factor is always positive, the remaining factor provides the condition for the monotonic increase:

$$ \frac{2}{\alpha} + 2\mu_0 > h. $$

For $\mu_0 > 0$, this condition is fulfilled because

$$ \frac{2}{\alpha} + 2\mu_0 > \mu_0 > h. $$

Here we used $\alpha > 0$, $\mu_0 > 0$, and $\mu_0 > h_0$. With $F(h) > 0$ strictly monotonically increasing and the linear function in

Eq. (11) monotonically decreasing there must be exactly one fixed-point solution $h_0 \in ] -\infty, \mu_0[$ with a positive firing rate $r_0 > 0$.

## APPENDIX J: SUMMARY TABLES

An overview of the parameters and the population models in different parameter regimes are provided in Tables II and III, respectively.

[1] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, Stochastic gene expression in a single cell, Science 297, 1183 (2002).

[2] A. Raj and A. Van Oudenaarden, Nature, nurture, or chance: Stochastic gene expression and its consequences, Cell 135, 216 (2008).

[3] S. Taheri-Araghi, S. Bradde, J. T. Sauls, N. S. Hill, P. A. Levin, J. Paulsson, M. Vergassola, and S. Jun, Cell-Size control and homeostasis in bacteria, Curr. Biol. 25, 385 (2015).

[4] S. Arzash, I. Tah, A. J. Liu, and M. L. Manning, Rigidity of epithelial tissues as a double optimization problem, Phys. Rev. Res. 7, 013157 (2025).

[5] A. A. Faisal, L. P. J. Selen, and D. M. Wolpert, Noise in the nervous system, Nat. Rev. Neurosci. 9, 292 (2008).

[6] B. Doiron, A. Litwin-Kumar, R. Rosenbaum, G. K. Ocker, and K. Josić, The mechanics of state dependent neural correlations, Nat. Neurosci. 19, 383 (2016).

[7] L. Waschke, N. A. Kloosterman, J. Obleser, and D. D. Garrett, Behavior needs neural variability, Neuron 109, 751 (2021).

[8] T. B. Kepler and T. C. Elston, Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representation, Biophys. J. 81, 3116 (2001).

[9] T. Guérin, J. Prost, and J.-F. Joanny, Bidirectional motion of motor assemblies and the weak-noise escape problem, Phys. Rev. E 84, 041901 (2011).

[10] C. Touya, T. Schwalger, and B. Lindner, Relation between cooperative molecular motors and active Brownian particles, Phys. Rev. E 83, 051913 (2011).

[11] R. Ma, G. S. Klindt, I. H. Riedel-Kruse, F. Jülicher, and B. M. Friedrich, Active phase and amplitude fluctuations of flagellar beating, Phys. Rev. Lett. 113, 048101 (2014).

[12] A. Dobrinevski and E. Frey, Extinction in neutrally stable stochastic Lotka-Volterra models, Phys. Rev. E 85, 051903 (2012).

[13] H. R. Wilson and J. D. Cowan, Excitatory and inhibitory interactions in localized populations of model neurons, Biophys. J. 12, 1 (1972).

[14] Y. Lu and J. Rinzel, Firing rate models for gamma oscillations in I-I and E-I networks, J. Comput. Neurosci. 52, 247 (2024).

[15] D. Levenstein, G. Buzsáki, and J. Rinzel, Nrem sleep in the rodent neocortex and hippocampus reflects excitable dynamics, Nat. Commun. 10, 2478 (2019).

[16] Y. Ahmadian, D. B. Rubin, and K. D. Miller, Analysis of the stabilized supralinear network, Neural Comput. 25, 1994 (2013).

[17] K.-F. Wong and X.-J. Wang, A recurrent network mechanism of time integration in perceptual decisions, J. Neurosci. 26, 1314 (2006).

[18] A. Shpiro, R. Moreno-Bote, N. Rubin, and J. Rinzel, Balance between noise and adaptation in competition models of perceptual bistability, J. Comput. Neurosci. 27, 37 (2009).

[19] G. Dumont, G. Northoff, and A. Longtin, Linear noise approximation for oscillations in a stochastic inhibitory network with delay, Phys. Rev. E 90, 012702 (2014).

[20] T. Kanashiro, G. K. Ocker, M. R. Cohen, and B. Doiron, Attentional modulation of neuronal variability in circuit models of cortex, eLife 6, e23978 (2017).

[21] G. Hennequin, Y. Ahmadian, D. B. Rubin, M. Lengyel, and K. D. Miller, The dynamical regime of sensory cortex: Stable dynamics around a single stimulus-tuned attractor account for patterns of noise variability, Neuron 98, 846 (2018).

[22] S.-I. Amari, Characteristics of random nets of analog neuron-like elements, IEEE Trans. Syst. Man Cybern. SMC-2, 643 (1972).

[23] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, Neural Dynamics: From Single Neurons to Networks and Models of Cognition (Cambridge University Press, Cambridge, UK, 2014).

[24] G. Deco, V. K. Jirsa, P. A. Robinson, M. Breakspear, and K. Friston, The dynamic brain: From spiking neurons to neural masses and cortical fields, PLoS Comput. Biol. 4, e1000092 (2008).

[25] O. D. Faugeras, J. D. Touboul, and B. Cessac, A constructive mean-field analysis of multi population neural networks with random synaptic weights and stochastic inputs, Front. Comput. Neurosci. 3, 323 (2009).

[26] S. El Boustani and A. Destexhe, A master equation formalism for macroscopic modeling of asynchronous irregular activity states, Neural Comput. 21, 46 (2009).

[27] G. La Camera, The mean field approach for populations of spiking neurons, in Computational Modelling of the Brain: Modelling Approaches to Cells, Circuits and Networks, edited by M. Giugliano, M. Negrello, and D. Lino (Springer International, Cham, 2022), pp. 125–157.

[28] G. K. Ocker, Republished: Dynamics of stochastic integrate-and-fire networks, Phys. Rev. X 13, 041047 (2023).

[29] N. Brunel and V. Hakim, Fast global oscillations in networks of integrate-and-fire neurons with low firing rates, Neural Comput. 11, 1621 (1999).

[30] N. Brunel, Sparsely connected networks of spiking neurons, J. Comput. Neurosci. 8, 183 (2000).
[31] W. Gerstner, Population dynamics of spiking neurons: Fast transients, asynchronous states, and locking, Neural Comput. **12**, 43 (2000).

[32] E. Montbrió, D. Pazó, and A. Roxin, Macroscopic description for networks of spiking neurons, Phys. Rev. X **5**, 021028 (2015).

[33] V. Schmutz, E. Löcherbach, and T. Schwalger, On a finite-size neuronal population equation, SIAM J. Appl. Dyn. Syst. **22**, 996 (2023).

[34] D. Q. Nykamp and D. Tranchina, A population density approach that facilitates large-scale modeling of neural networks: Analysis and an application to orientation tuning, J. Comput. Neurosci. **8**, 19 (2000).

[35] N. Cain, R. Iyer, C. Koch, and S. Mihalas, The computational properties of a simplified cortical column model, PLoS Comput. Biol. **12**, e1005045 (2016).

[36] T. Schwalger, M. Deger, and W. Gerstner, Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size, PLoS Comput. Biol. **13**, e1005507 (2017).

[37] M. Augustin, J. Ladenbauer, F. Baumann, and K. Obermayer, Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: Comparison and implementation, PLoS Comput. Biol. **13**, e1005545 (2017).

[38] A. René, A. Longtin, and J. H. Macke, Inference of a mesoscopic population model from population spike trains, Neural Comput. **32**, 1448 (2020).

[39] S. Wang, V. Schmutz, G. Bellec, and W. Gerstner, Mesoscopic modeling of hidden spiking neurons, Adv. Neural Inf. Process. Syst. **35**, 23566 (2022).

[40] C. Pozzorini, S. Mensi, O. Hagens, R. Naud, C. Koch, and W. Gerstner, Automated high-throughput characterization of single neurons by means of simplified spiking models, PLoS Comput. Biol. **11**, e1004275 (2015).

[41] Y. N. Billeh, B. Cai, S. L. Gratiy, K. Dai, R. Iyer, N. W. Gouwens, R. Abbasi-Asl, X. Jia, J. H. Siegle, S. R. Olsen, C. Koch, S. Mihalas, and A. Arkhipov, Systematic integration of structural and functional data into multi-scale models of mouse primary visual cortex, Neuron **106**, 388 (2020).

[42] M. A. Buice, J. D. Cowan, and C. C. Chow, Systematic fluctuation expansion for neural network activity equations, Neural Comput. **22**, 377 (2010).

[43] M. A. Buice and C. C. Chow, Dynamic finite size effects in spiking neural networks, PLoS Comput. Biol. **9**, e1002872 (2013).

[44] G. Dumont, A. Payeur, and A. Longtin, A stochastic-field description of finite-size spiking neural networks, PLoS Comput. Biol. **13**, e1005691 (2017).

[45] S. Heesen and W. Stannat, Fluctuation limits for mean-field interacting nonlinear Hawkes processes, Stochast. Process. Appl. **139**, 280 (2021).

[46] V. V. Klinshov and S. Y. Kirillov, Shot noise in next-generation neural mass models for finite-size networks, Phys. Rev. E **106**, L062302 (2022).

[47] M. Deger, T. Schwalger, R. Naud, and W. Gerstner, Fluctuations and information filtering in coupled populations of spiking neurons with adaptation, Phys. Rev. E **90**, 062704 (2014).

[48] H. Sompolinsky, A. Crisanti, and H. J. Sommers, Chaos in random neural networks, Phys. Rev. Lett. **61**, 259 (1988).

[49] L. Molgedey, J. Schuchhardt, and H. G. Schuster, Suppressing chaos in neural networks by noise, Phys. Rev. Lett. **69**, 3717 (1992).

[50] C. van Vreeswijk and H. Sompolinsky, Chaos in neuronal networks with balanced excitatory and inhibitory activity, Science **274**, 1724 (1996).

[51] A. Renart, R. Moreno-Bote, X.-J. Wang, and N. Parga, Mean-driven and fluctuation-driven persistent activity in recurrent networks, Neural Comput. **19**, 1 (2007).

[52] M. Helias, T. Tetzlaff, and M. Diesmann, The correlation structure of local neuronal networks intrinsically results from recurrent dynamics, PLoS Comput. Biol. **10**, e1003428 (2014).

[53] D. S. Goldobin, M. di Volo, and A. Torcini, Reduction methodology for fluctuation driven population dynamics, Phys. Rev. Lett. **127**, 038301 (2021).

[54] S. Lefort, C. Tomm, J. C. F. Sarria, and C. C. H. Petersen, The excitatory neuronal network of the c2 barrel column in mouse primary somatosensory cortex, Neuron **61**, 301 (2009).

[55] A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang, Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model, Cereb. Cortex **10**, 910 (2000).

[56] L. Mazzucato, A. Fontanini, and G. La Camera, Dynamics of multistable states during ongoing and evoked cortical activity, J. Neurosci. **35**, 8214 (2015).

[57] G. Gigante, N. Deco, S. Marom, and P. Del Giudice, Network events on multiple space and time scales in cultured neural networks and in a stochastic rate model, PLoS Comput. Biol. **11**, e1004547 (2015).

[58] V. Schmutz, W. Gerstner, and T. Schwalger, Mesoscopic population equations for spiking neural networks with synaptic short-term plasticity, J. Math. Neurosc. **10**, 5 (2020).

[59] 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**, e1010809 (2022).

[60] H. Ko, S. B. Hofer, B. Pichler, K. A. Buchanan, P. J. Sjöström, and T. D. Mrsic-Flogel, Functional specificity of local synaptic connections in neocortical networks, Nature (London) **473**, 87 (2011).

[61] M. M. Churchland, M. Y. Byron, J. P. Cunningham, L. P. Sugrue, M. R. Cohen, G. S. Corrado, W. T. Newsome, A. M. Clark, P. Hosseini, B. B. Scott et al., Stimulus onset quenches neural variability: A widespread cortical phenomenon, Nat. Neurosci. **13**, 369 (2010).

[62] S. Wieland, D. Bernardi, T. Schwalger, and B. Lindner, Slow fluctuations in recurrent networks of spiking neurons, Phys. Rev. E **92**, 040901(R) (2015).

[63] A. Sanzeni, M. H. Histed, and N. Brunel, Response nonlinearities in networks of spiking neurons, PLoS Comput. Biol. **16**, e1008165 (2020).

[64] N. Schieferstein, T. Schwalger, B. Lindner, and R. Kempter, Intra-ripple frequency accommodation in an inhibitory network model for hippocampal ripple oscillations, PLoS Comput. Biol. **20**, e1011886 (2024).

[65] P. Ekelmans, N. Kraynyukova, and T. Tchumatchenko, Targeting operational regimes of interest in recurrent neural networks, PLoS Comput. Biol. **19**, e1011097 (2023).

[66] M. Mattia and P. Del Giudice, Population dynamics of interacting spiking neurons, Phys. Rev. E **66**, 051917 (2002).
[67] B. Derrida and Y. Pomeau, Random networks of automata: A simple annealed approximation, *Europhys. Lett.* **1**, 45 (1986).

[68] R. Jolivet, A. Rauch, H. R. Lüscher, and W. Gerstner, Predicting spike timing of neocortical pyramidal neurons by simple threshold models, *J. Comput. Neurosci.* **21**, 35 (2006).

[69] N. J. Priebe and D. Ferster, Inhibition, spike threshold, and stimulus selectivity in primary visual cortex, *Neuron* **57**, 482 (2008).

[70] N. Brunel and D. Hansel, How noise affects the synchronization properties of recurrent networks of inhibitory neurons, *Neural Comput.* **18**, 1066 (2006).

[71] B. Lindner, B. Doiron, and A. Longtin, Theory of oscillatory firing induced by spatially correlated noise and delayed inhibitory feedback, *Phys. Rev. E* **72**, 061919 (2005).

[72] S. P. Muscinelli, W. Gerstner, and T. Schwalger, How single neuron properties shape chaotic dynamics and signal transmission in random neural networks, *PLoS Comput. Biol.* **15**, e1007122 (2019).

[73] N. Brunel and V. Hakim, Sparsely synchronized neuronal oscillations, *Chaos* **18**, 015113 (2008).

[74] E. Ledoux and N. Brunel, Dynamics of networks of excitatory and inhibitory neurons in response to time-dependent inputs, *Front. Comput. Neurosci.* **5**, 25 (2011).

[75] J. Kadmon and H. Sompolinsky, Transition to chaos in random neuronal networks, *Phys. Rev. X* **5**, 041030 (2015).

[76] D. Dahmen, S. Grün, M. Diesmann, and M. Helias, Second type of criticality in the brain uncovers rich multiple-neuron dynamics, *Proc. Natl. Acad. Sci. USA* **116**, 13051 (2019).

[77] G. Hennequin, T. P. Vogels, and W. Gerstner, Optimal control of transient dynamics in balanced networks supports generation of complex movements, *Neuron* **82**, 1394 (2014).

[78] Y.-L. Shi, R. Zeraati, I. B. Laboratory, A. Levina, and T. A. Engel, Brain-wide organization of intrinsic timescales at single-neuron resolution, https://doi.org/10.1101/2025.08.30.673281 (2025).

[79] A. Kulkarni, J. Ranft, and V. Hakim, Synchronization, stochasticity, and phase waves in neuronal networks with spatially-structured connectivity, *Front. Comput. Neurosci.* **14**, 569644 (2020).

[80] A. Sanzeni, B. Akitake, H. C. Goldbach, C. E. Leedy, N. Brunel, and M. H. Histed, Inhibition stabilization is a widespread property of cortical networks, *eLife* **9**, e54875 (2020).

[81] L. Kang, J. Ranft, and V. Hakim, Beta oscillations and waves in motor cortex can be accounted for by the interplay of spatially structured connectivity and fluctuating inputs, *eLife* **12**, e81446 (2023).

[82] G. Hahn, A. Kumar, H. Schmidt, T. R. Knösche, and G. Deco, Rate and oscillatory switching dynamics of a multilayer visual microcircuit model, *eLife* **11**, e77594 (2022).

[83] L. Yu and T. O. Taillefumier, Metastable spiking networks in the replica-mean-field limit, *PLoS Comput. Biol.* **18**, e1010215 (2022).

[84] G. Dumont, J. Henry, and C. O. Tarniceriu, Oscillations in a fully connected network of leaky integrate-and-fire neurons with a poisson spiking mechanism, *J. Nonlin. Sci.* **34**, 18 (2024).

[85] B. W. Knight, D. Manin, and L. Sirovich, Dynamical models of interacting neuron populations in visual cortex, in *Symposium on Robotics and Cybernetics: Computational Engineering in Systems Applications*, edited by E. C. Geil (Cite Scientifique, Lille, 1996), Vol. 54, pp. 4–8.

[86] G. V. Vinci, R. Benzi, and M. Mattia, Self-consistent stochastic dynamics for finite-size networks of spiking neurons, *Phys. Rev. Lett.* **130**, 097402 (2023).

[87] https://github.com/Nils-Greven/random-connectivity-finite-size-network.

[88] H. Risken, *The Fokker-Planck Equation* (Springer, Berlin, 1984).

[89] D. B. Owen, A table of normal integrals, *Commun. Stat. Simul. Comput.* **9**, 389 (1980).