m o#2#1#1_#2 m o#2C_#1C_#1^#2 m o#2˙#1˙#1_#2 m o#2˙C_#1˙C_#1^#2 m oE_#1^#2 m
o˙E_#1^#2

# The determining role of covariances in large networks of stochastic neurons

Vincent Painchaud  Department of Mathematics and Statistics, McGill
University, Montreal, Québec, Canada H3A 0B6  Patrick Desrosiers  Nicolas
Doyon

Biological neural networks are notoriously hard to model due to their
stochastic behavior and high dimensionality. We tackle this problem by
constructing a dynamical model of both the expectations and covariances of the
fractions of active and refractory neurons in the network’s populations. We do
so by describing the evolution of the states of individual neurons with a
continuous-time Markov chain, from which we formally derive a low-dimensional
dynamical system. This is done by solving a moment closure problem in a way
that is compatible with the nonlinearity and boundedness of the activation
function. Our dynamical system captures the behavior of the high-dimensional
stochastic model even in cases where the mean-field approximation fails to do
so. Taking into account the second-order moments modifies the solutions that
would be obtained with the mean-field approximation, and can lead to the
appearance or disappearance of fixed points and limit cycles. We moreover
perform numerical experiments where the mean-field approximation leads to
periodically oscillating solutions, while the solutions of the second-order
model can be interpreted as an average taken over many realizations of the
stochastic model. Altogether, our results highlight the importance of
including higher moments when studying stochastic networks and deepen our
understanding of correlated neuronal activity.

## 1 Introduction

Neuronal activity is intrinsically stochastic. Sources of randomness have been
identified at different scales, ranging from spontaneous neurotransmitters
release to noise amplification at the network level, and have been shown to
contribute to cellular and behavioral trial-to-trial variability [13, 56, 22].
Although stochasticity has often been treated as a nuisance, recent studies
have highlighted its determining role in neural population coding [2] and in
many brain functions, such as awareness [60], decision making [37], and
working memory [52]. Moreover, biologically realistic numerical simulations
have demonstrated that stochastic networks can support spike-time coding with
millisecond precision [43] while abstract models have revealed the
computational advantage of stochastic spiking neurons over deterministic ones
[41].

From a mathematical point of view, incorporating stochasticity into neural
network models gives rise to major difficulties, especially when trying to
make interpretable predictions about their large-scale dynamical behavior. To
address this, a common practice is to perform mean-field approximations [8,
40, 21, 31], which have long been used in statistical mechanics when studying
phase transitions in spin systems [19, 38]. Perhaps the most important example
of a mean-field model is the one developed by Wilson and Cowan in 1972 [58],
building on the work of Beurle [3] and Griffith [29, 30]. This model achieved
immediate success due to its ability to represent hysteresis phenomena as well
as oscillations in biological networks of excitatory and inhibitory neurons.
It is still widely used today, sometimes as a firing-rate model [21, 26, 55,
33], and has been the starting point of several extensions [18, 7, 17, 15,
57].

As it is based on a mean-field approximation, the Wilson–Cowan model is
restricted to modeling only the average behavior of neural networks, and
cannot represent correlations between the states of different neurons.
However, correlations have been found to be important in brain activity [50,
53, 1, 49, 2], and there has been ongoing effort to generalize the
Wilson–Cowan model to take them into account [15]. An important example of
such work is that started by Buice and Cowan [10, 11], later joined by Chow
[12, 9] and followed by Bressloff [6, 5], who used path integral methods to
derive a master equation that generalizes the Wilson–Cowan model. This
approach has been used to study the relationship between structure and
activity in neural networks [45].

Here, we propose a new way to generalize the Wilson–Cowan model to include
covariances between the activities of neural populations. The idea is to first
build a stochastic process to model the dynamics from the microscopic point of
view of individual neurons, and then to use this process to obtain a
nonautonomous differential equation that models the same dynamics, but from
the macroscopic point of view of neural populations. This nonautonomous
differential equation poses a moment closure problem (see [36] for an
extensive discussion), which can be solved while still considering both first-
and second-order moments.

The moment closure framework has not been very popular in computational
neuroscience so far, but it has been extensively used in other areas of
mathematical biology during the last decades. In particular, it has been used
to include covariances in compartmental models in epidemiology [39, 34, 32,
54, 14, 35] as well as in population dynamics models in ecology [42, 51]. One
of the main advantages of this method is that it reveals precisely where
covariances can have an effect on the dynamics. Moreover, it allows a more
systematic treatment of the neurons’ refractory period, as it has been done in
[48].

The paper is organized as follows. First, we construct a continuous-time
Markov chain that describes the evolution of the state of each node in a large
network, the goal being to mimic the behavior of biological neurons. This
Markov chain is similar to one already proposed (but not extensively studied)
by Cowan [16] and is reminiscent of a process recently studied by Zarepour et
al. [59]. Then, we split the network into a small set of large populations and
derive a nonautonomous system of differential equations that describes the
evolution of the expected state of each population. Using a new approximation
of the expectation of a sigmoid function, we find a solution to the moment
closure problem that involves covariances between population state variables.
Finally, we present three examples that demonstrate the impacts of
covariances. The first one illustrates how covariances can increase the
accuracy when predicting the macroscopic behavior of the Markov chain. In the
other two examples, the impact of including covariances in the dynamical
system goes beyond error reduction: the second-order system reaches a steady
state where covariances are nonzero, implying that this steady state is
inaccessible to the mean-field model.

## 2 The model

We seek to describe the dynamics of a large biological neural network from a
macroscopic point of view in such a way that the resulting model generalizes
Wilson–Cowan’s, but also includes correlations between states of neurons. We
are also looking for a description that would depend on statistical properties
of neuron parameters rather than on their precise values. To do so, we
construct a Markov chain to describe the states of all neurons of the network.
Then, we split the network into populations and we derive a dynamical system
that describes the evolution of the states of populations. The model
constructed in this section is also presented (under slightly different
assumptions) in [48], and it is a special case of a more general model
presented in [46, Chapter 2].

### 2.1 A modelization of a biological neural network

We consider a network of $N$ neurons. Links between neurons are described by a
random real-valued $N\times N$ matrix $W$ that we call the _weight matrix_ ,
and that is defined on a probability space $(H,\mathscr{H},\mu)$. An entry
$W_{jk}$ of the weight matrix describes the connection from neuron $k$ to
neuron $j$; it can be either positive or negative, describing either an
excitatory or an inhibitory connection.

To study the network’s dynamics from a macroscopic point of view, we suppose
that the neurons are split into a small number $n$ of distinct subgroups, each
consisting of a large number of neurons, which we call _populations_. This
splitting is made through a partition $\mathscr{P}$ of the set
$\\{1,\ldots,N\\}$ where each $J\in\mathscr{P}$ corresponds to a population.
We assume that the weights $W_{jk}$ are independent random variables,
identically distributed over populations, and assume the same for all
parameters that will be introduced in the next section.

### 2.2 A description of the evolution of the network’s state

We want the dynamics of our network to model the behavior of biological
neurons. The basic behavior we seek to represent is that a neuron fires when
it receives sufficient input, and that after firing occurs a short refractory
period during which a neuron cannot fire again. To mimic this, we assume that
the nodes of our network can take the states described on Fig. 1, where the
active state corresponds to that of a neuron which is firing. The whole
network’s state can then be described by an element of the set
$E\vcentcolon=\\{0,1,i\\}^{N}$.

$0$$1$$i$$\alpha$$\beta$$\gamma$

$0$ | sensitive state
---|---
$1$ | active state
$i$ | refractory state
$\alpha$ | activation rate with sufficient input
$\beta$ | deactivation rate
$\gamma$ | sensitization rate

Figure 1: Neurons’ states and allowed transitions between them with
corresponding rates. Here, $i$ denotes the imaginary unit.

To describe the transitions between the states, we consider for each neuron
$j$ three positive random variables $\alpha_{j}$, $\beta_{j}$ and $\gamma_{j}$
and a real-valued one $\theta_{j}$, which are all independent and defined on
$(H,\mathscr{H},\mu)$. The possible transitions are described on Fig. 1. The
parameters $\beta_{j}$ and $\gamma_{j}$ are the rates at which neuron $j$ can
make the transitions $1\mapsto i$ and $i\mapsto 0$. The activation rate is
described by a function $a_{j}\colon H\times E\to(0,\infty)$ given by

$a_{j}(\eta,x)\vcentcolon=\alpha_{j}(\eta){}_{T_{j}(x)}(\eta),$

where ${}_{T_{j}(x)}$ is the indicator function of the set

$T_{j}(x)\vcentcolon=\Bigl{\\{}\eta\in
H:\sum_{k=1}^{N}W_{jk}(\eta)\operatorname{Re}x_{k}+Q_{J}>\theta_{j}(\eta)\Bigr{\\}},$

$Q_{J}\in\mathbb{R}$ being a deterministic external input received by the
population $J\in\mathscr{P}$ to which belongs neuron $j$. Thus, $j$ can
activate at rate $\alpha_{j}$ when its input exceeds its threshold
$\theta_{j}$.

This description of the microscopic behavior is the intuitive description of a
continuous-time Markov chain, which we now define properly. We introduce for
each choice of parameters $\eta\in H$ a matrix
$M^{\eta}=\\{m^{\eta}(x,y):x,y\in E\\}$ with entries

$m^{\eta}(x,y)\vcentcolon=\sum_{j=1}^{N}m_{j}^{\eta}(x,y)\smash{\prod_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{N}}\delta_{x_{k}y_{k}}$

where

$\displaystyle m_{j}^{\eta}(x,y)\vcentcolon=a_{j}(\eta,x)$
$\displaystyle(1-\lvert x_{j}\rvert)\bigl{(}\operatorname{Re}y_{j}-(1-\lvert
y_{j}\rvert)\bigr{)}$ $\displaystyle{}+\beta_{j}(\eta)$
$\displaystyle\operatorname{Re}(x_{j})(\operatorname{Im}y_{j}-\operatorname{Re}y_{j})$
$\displaystyle{}+\gamma_{j}(\eta)$
$\displaystyle\operatorname{Im}(x_{j})\bigl{(}(1-\lvert
y_{j}\rvert)-\operatorname{Im}y_{j}\bigr{)}.$

To make sense of $m^{\eta}_{j}(x,y)$, recall that $x_{j}\in\\{0,1,i\\}$, so
that exactly one of $\operatorname{Re}x_{j}$, $\operatorname{Im}x_{j}$ and
$1-\lvert x_{j}\rvert$ is $1$ while the other two are $0$, so that
$m^{\eta}_{j}(x,y)$ is always $\pm a_{j}(\eta,x)$, $\pm\beta_{j}(\eta)$ or
$\pm\gamma_{j}(\eta)$. A simple calculation then shows that the matrix
$M^{\eta}$ is the generator of a continuous-time Markov chain. Thus, it
follows from the Kolmogorov extension theorem (see e.g. [20, 44] for details)
that there exists a probability measure $\mathbb{P}^{\eta}$ on
$(\Omega,\mathscr{F})\vcentcolon=(E,2^{E})^{[0,\infty)}$ such that for any
$x,y\in E$, as $\Delta t\downarrow 0$,

$\mathbb{P}^{\eta}[X_{t+\Delta
t}=y\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]=\delta_{xy}+m^{\eta}(x,y)\Delta
t+o(\Delta t),$

where $\\{X_{t}\\}_{t\geq 0}$ is the coordinate mapping process
$X_{t}(\omega)\vcentcolon=\omega(t)$ on $(\Omega,\mathscr{F})$. In particular,
if $x\in E$ has $x_{j}=0$, then

$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]$
$\displaystyle=a_{j}(\eta,x)\Delta t+o(\Delta t),$ (1a) while in general
$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}^{j}=1]$
$\displaystyle=\beta_{j}(\eta)\Delta t+o(\Delta t),$ (1b)
$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=0\nonscript\,|\nonscript\,\mathopen{}X_{t}^{j}=i]$
$\displaystyle=\gamma_{j}(\eta)\Delta t+o(\Delta t),$ (1c)

and the forbidden transitions $0\mapsto i$, $i\mapsto 1$ and $1\mapsto 0$ all
have $o(\Delta t)$ rates. Hence, the stochastic process $\\{X_{t}\\}_{t\geq
0}$ describes the state of the network with respect to time as described at
the beginning of the section.

### 2.3 A macroscopic approximation of the Markov chain’s behavior

The above Markov chain completely describes the dynamics of the network from a
microscopic point of view. Now, we want to find an approximation of this
Markov chain to describe the macroscopic behavior of the network. To make this
goal more precise, we introduce for each population $J\in\mathscr{P}$ the
processes

$\displaystyle A_{t}^{J}$ $\displaystyle\vcentcolon=\frac{1}{\lvert
J\rvert}\sum_{j\in J}\operatorname{Re}X_{t}^{j},$ (2a) $\displaystyle
R_{t}^{J}$ $\displaystyle\vcentcolon=\frac{1}{\lvert J\rvert}\sum_{j\in
J}\operatorname{Im}X_{t}^{j},$ (2b) $\displaystyle S_{t}^{J}$
$\displaystyle\vcentcolon=\frac{1}{\lvert J\rvert}\sum_{j\in
J}\bigl{(}1-\lvert X_{t}^{j}\rvert\bigr{)},$ (2c)

which are the fractions of neurons in population $J$ that are in the active,
refractory and sensitive state respectively. Our goal is to find a dynamical
system to describe the evolution of

$\E{A}[J](t)\vcentcolon=\mathbb{E}^{\eta}[A_{t}^{J}],\quad\E{R}[J](t)\vcentcolon=\mathbb{E}^{\eta}[R_{t}^{J}],\quad\E{S}[J](t)\vcentcolon=\mathbb{E}^{\eta}[S_{t}^{J}],$
(3)

where $\mathbb{E}^{\eta}$ denotes the expectation on
$(\Omega,\mathscr{F},\mathbb{P}^{\eta})$.

A first step towards an understanding of the evolution of these expected
fractions of populations is a description of the evolution of the
probabilities

$\displaystyle p_{j}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=1]$
$\displaystyle=\mathbb{E}^{\eta}\big{[}\operatorname{Re}X_{t}^{j}\big{]},$
$\displaystyle r_{j}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=i]$
$\displaystyle=\mathbb{E}^{\eta}\big{[}\operatorname{Im}X_{t}^{j}\big{]},$
$\displaystyle q_{j}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=0]$
$\displaystyle=\mathbb{E}^{\eta}\big{[}1-\lvert X_{t}^{j}\rvert\big{]}.$

By considering these probabilities at a time $t+\Delta t$ and conditioning
over all events $\\{X_{t}=x\\}$ for $x\in E$, it is not hard to verify using
the transition probabilities given in eq. (1) that

$\displaystyle\dot{p}_{j}(t)$
$\displaystyle=-\beta_{j}(\eta)p_{j}(t)+\mathbb{E}^{\eta}\big{[}a_{j}(\eta,X_{t})\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\big{]},$ (4a) $\displaystyle\dot{q}_{j}(t)$
$\displaystyle=-\mathbb{E}^{\eta}\big{[}a_{j}(\eta,X_{t})\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\big{]}+\gamma_{j}(\eta)r_{j}(t),$ (4b)
$\displaystyle\dot{r}_{j}(t)$
$\displaystyle=-\gamma_{j}(\eta)r_{j}(t)+\beta_{j}(\eta)p_{j}(t).$ (4c)

The details are given in Appendix A.1. We use these differential equations as
the starting point to describe the evolution of the expected fractions of
populations defined in eq. (3).

To obtain expressions for the derivatives of $\E{A}[J]$, $\E{R}[J]$ and
$\E{S}[J]$, the idea is to average the derivatives given in eq. (4) over $J$
using the linearity of expectation and derivatives. For example,

$\dE{A}[J](t)=\mathbb{E}^{\eta}\bigg{[}\frac{1}{\lvert J\rvert}\sum_{j\in
J}\Bigl{(}-\beta_{j}(\eta)\operatorname{Re}X_{t}^{j}+a_{j}(\eta,X_{t})\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\Bigr{)}\bigg{]}.$

Assuming that each population is large, the law of large numbers now motivates
to approximate the averages of transition rates over populations by their
expectation. Thus, $\frac{1}{\lvert J\rvert}\sum_{j\in
J}\beta_{j}(\eta)\operatorname{Re}X_{t}^{j}\approx\beta_{J}A_{t}^{J}$ where
$\beta_{J}\vcentcolon=\mathbb{E}_{\mu}[\beta_{j}]$ for $j\in J$,
$\mathbb{E}_{\mu}$ being the expectation on $(H,\mathscr{H},\mu)$. The
activation rates are handled in the same way, but since these rates are step
functions, their expectation results in the cumulative distribution function
$F_{\theta_{J}}$ of the thresholds in population $J$ evaluated at a
population-averaged input:

$\frac{1}{\lvert J\rvert}\sum_{j\in J}a_{j}(\eta,X_{t})\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\approx\alpha_{J}F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}$

where $\alpha_{J}\vcentcolon=\mathbb{E}_{\mu}[\alpha_{j}]$ and

$B_{t}^{J}\vcentcolon=\sum_{K\in\mathscr{P}}c_{JK}A_{t}^{K}+Q_{J}\quad\text{with}\quad
c_{JK}\vcentcolon=\lvert K\rvert\mathbb{E}_{\mu}[W_{jk}]$

for $j\in J$ and $k\in K$. The complete details, including the limits of these
approximations, are given in Appendix A.

This method results in approximate expressions for the derivatives of
$\E{A}[J]$, $\E{R}[J]$ and $\E{S}[J]$, and ultimately leads to model the
macroscopic dynamics of the network by the differential equations

$\displaystyle\dE{A}[J](t)$
$\displaystyle=-\beta_{J}\E{A}[J](t)+\alpha_{J}\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}],$
(5a) $\displaystyle\dE{R}[J](t)$
$\displaystyle=-\gamma_{J}\E{R}[J](t)+\beta_{J}\E{A}[J](t),$ (5b)
$\displaystyle\dE{S}[J](t)$
$\displaystyle=-\alpha_{J}\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]+\gamma_{J}\E{R}[J](t),$
(5c)

where $\gamma_{J}\vcentcolon=\mathbb{E}_{\mu}[\gamma_{j}]$ for $j\in J$. For
each population, one of these three equations is redundant since
$A_{t}^{J}+R_{t}^{J}+S_{t}^{J}\equiv 1$ for all $t$. Thus, if the network has
$n$ populations, eq. (5) corresponds to $2n$ independent differential
equations. In the following, we use the active and refractory fractions as the
independent variables, and always see sensitive fractions simply as a function
of them.

A crucial aspect of eq. (5) is that it is not autonomous: there is an explicit
time dependence in the expectation
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$, which is additionally
an unknown function. Further assumptions or approximations are required to
obtain an autonomous system—this corresponds to a moment closure problem [36].
The simplest way to close the system would be to use the mean-field
approximation $\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]\approx
F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}\E{S}[J](t)$, where
$\E{B}[J](t)\vcentcolon=\mathbb{E}^{\eta}[B_{t}^{J}]$. This results in the
dynamical system

$\displaystyle\dE{A}[J]$
$\displaystyle=-\beta_{J}\E{A}[J]+\alpha_{J}F_{\theta_{J}}(\E{B}[J])\E{S}[J],$
(6a) $\displaystyle\dE{R}[J]$
$\displaystyle=-\gamma_{J}\E{R}[J]+\beta_{J}\E{A}[J],$ (6b)

which is closed since for each population $J$, $\E{S}[J]$ and $\E{B}[J]$ are
functions of expectations of active and refractory fractions. Since the
approximation used to close the system does not depend on the specific choice
of parameters $\eta\in H$, any explicit dependence on $\eta$ disappeared.
Thus, the macroscopic dynamics predicted by the model only depend on the
averages of parameters over populations. This will be the case for any
macroscopic model obtained from a closure of the system given in eq. (5), as
long as the approximations used to close the system do not depend explicitly
on $\eta$.

The mean-field system given in eq. (6) is studied in details in [48] where it
is shown to be closely related to Wilson–Cowan’s equations [58]: if the
refractory fractions are forced to their equilibrium solutions
$\E{R}[J]=\frac{\beta_{J}}{\gamma_{J}}\E{A}[J]$, the resulting system is
equivalent to Wilson–Cowan’s. In the next section, we go further and present a
solution to the moment closure problem posed by eq. (5) that includes second-
order moments.

## 3 A second-order solution to the moment closure problem

We are looking for a solution to the moment closure problem posed by eq. (5)
that includes covariances between fractions of populations. The first step is
to find how the second moments of the active and refractory fractions of
populations evolve in time. To simplify notation in the following, we define

$\C{YZ}[JK](t)\vcentcolon=\mathrm{Cov}[Y_{t}^{J},Z_{t}^{K}]=\mathbb{E}^{\eta}[Y_{t}^{J}Z_{t}^{k}]-\mathbb{E}^{\eta}[Y_{t}^{J}]\mathbb{E}^{\eta}[Z_{t}^{K}],$

where $Y$ and $Z$ stand for either $A$, $R$, $S$ or $B$, and
$J,K\in\mathscr{P}$. Then, with the method we used above to find eq. (5), it
can be shown (and this is done in Appendix A) that

$\displaystyle\begin{split}\dC{AA}[JK](t)&=-(\beta_{J}+\beta_{K})\C{AA}[JK](t)\\\
&\qquad+\alpha_{K}\mathrm{Cov}[A_{t}^{J},F_{\theta_{K}}(B_{t}^{K})S_{t}^{K}]\\\
&\qquad+\alpha_{J}\mathrm{Cov}[A_{t}^{K},F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}],\end{split}$
(7a)
$\displaystyle\begin{split}\dC{RR}[JK](t)&=-(\gamma_{J}+\gamma_{K})\C{RR}[JK](t)\\\
&\qquad+\beta_{K}\C{AR}[KJ](t)+\beta_{J}\C{AR}[JK](t),\end{split}$ (7b)
$\displaystyle\begin{split}\dC{AR}[JK](t)&=-(\beta_{J}+\gamma_{K})\C{AR}[JK](t)\\\
&\qquad+\beta_{K}\C{AA}[JK](t)+\alpha_{J}\mathrm{Cov}[R_{t}^{K},F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}].\end{split}$
(7c)

Thus, to obtain a closed dynamical system involving only first and second
moments of active and refractory fractions, we need to find approximations to
expectations of the forms
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$,
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]$ and
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}R_{t}^{K}]$ in terms of
expectations and covariances of active and refractory fractions of
populations. Remark that, since sensitive fractions and inputs are expressed
in terms of active and refractory fractions by linear relations, the
bilinearity of covariance implies that any covariance of the form $\C{YZ}[JK]$
with $Y$ and $Z$ being $A$, $R$, $S$ or $B$ can be expressed using active and
refractory fractions only.

The simplest approach would be to assume that $F_{\theta_{J}}$ is smooth
enough to use a Taylor approximation, which would allow to express the
problematic expectations in terms of central moments of $A_{t}^{J}$ and
$R_{t}^{J}$. Then, neglecting all central moments of order 3 and higher, we
could get a closed dynamical system. This approach has already been taken in
[46], where the resulting dynamical system was studied in detail. In the case
of a single population, if the refractory fraction of the network is
neglected, the system also is similar to one that has been obtained in [12,
Eqs. (5.3) and (5.4)]. Unfortunately, the analysis of the resulting model [46,
Section 4.4] shows that it is not better than the mean-field model given by
eq. (6), because it does not predict more accurately the behavior of the
network, and because the covariances can blow up, in which case the model
loses all biological meaning. It is not hard to understand why this approach
fails: the expectations that need to be approximated are bounded functions,
whereas their approximations obtained from a Taylor expansion of
$F_{\theta_{J}}$ are not. Here, we propose another solution to this moment
closure problem built around a function that approximates the expectation of a
sigmoid function of a random variable in terms of its first and second
moments.

### 3.1 Approximation of a sigmoid function’s expectation

The difficulty in approximating
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$ or
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]$ comes, to a
large extent, from the cumulative distribution function $F_{\theta_{J}}$. As a
first step, we construct an approximation of the expectation of
$F_{\theta_{J}}(B_{t}^{J})$ only. To do so, we assume that in each population,
the thresholds follow a symmetric, unimodal distribution, so that
$F_{\theta_{J}}$ has a sigmoid shape with inflection point at the mean
$\theta_{J}$.

We now ask a simple question: if the random variable $B_{t}^{J}$ follows a
symmetric distribution with mean $\E{B}[J](t)$ and a small variance
$\C{BB}[JJ](t)$, how should the expectation of $F_{\theta_{J}}(B_{t}^{J})$
compare to $F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}$? It is not hard to find
a qualitative answer. Indeed, when $\E{B}[J](t)<\theta_{J}$, the fact that
$F_{\theta_{J}}$ increases faster to the right of $\E{B}[J](t)$ than it
decreases to its left has the effect that
$F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}$ must underestimate the actual
expectation of $F_{\theta_{J}}(B_{t}^{J})$. This effect is reversed when
$\E{B}[J](t)>\theta_{J}$. This is sketched on Fig. 2.

Figure 2: Example of comparison between $F_{\theta_{J}}(\E{B}[J](t))$, the
actual expectation $\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]$, and our
approximation $G_{J}(\E{B}[J](t),\C{BB}[JJ](t))$ for a fixed variance
$\C{BB}[JJ](t)$, assuming the thresholds follow a logistic distribution with
mean $\theta_{J}$ and scaling factor $s_{\theta_{J}}$. To compute
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]$, we assumed as a heuristic
device that $B_{t}^{J}$ follows a logistic distribution with mean
$\E{B}[J](t)$ and scaling factor $s_{\theta_{J}}$ as well, but this choice was
made for illustration purposes only.

Observe from Fig. 2 that the deformation of
$F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}$ into
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]$ resembles an increase in the
standard deviation of the distribution of which $F_{\theta_{J}}$ is the
cumulative distribution function. If $\sigma_{\theta_{J}}$ denotes the
standard deviation of a threshold $\theta_{j}$ of population $J$, then
increasing the standard deviation of $\theta_{j}$ by a factor $1+a$ for some
$a>0$ corresponds to changing $\theta_{j}$ to $\theta_{j}(1+a)-\theta_{J}a$.
The cumulative distribution functions of the two cases can be related as

$\mathbb{P}^{\eta}[\theta_{j}(1+a)-\theta_{J}a\leq
x]=\mathbb{P}^{\eta}\Big{[}\theta_{j}\leq\frac{x+\theta_{J}a}{1+a}\Big{]}.$

Therefore, our previous observations suggest to look for an approximation of
the form

$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]\approx
G_{J}\bigl{(}\E{B}[J](t),\C{BB}[JJ](t)\bigr{)}$ (8)

with

$G_{J}(b,v)\vcentcolon=F_{\theta_{J}}\Bigl{(}\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\Bigr{)},$
(9)

for some $g_{J}\colon\mathbb{R}\times[0,\infty)\to[0,\infty)$ to be
determined.

We can find a candidate for $g_{J}$ based on what $G_{J}$ should be for small
variances. Indeed, if the variance $\C{BB}[JJ](t)$ of $B_{t}^{J}$ is small,
then a second-order Taylor approximation of $F_{\theta_{J}}$ should provide a
reasonable approximation of $\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]$, so
that

$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]\approx
F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}+\frac{1}{2}F_{\theta_{J}}^{\prime\prime}\bigl{(}\E{B}[J](t)\bigr{)}\C{BB}[JJ](t).$
(10)

For the approximation given by eq. (8) to be consistent with this, it must be
that

$G_{J}(b,0)=F_{\theta_{J}}(b)\quad\text{and}\quad\partial_{2}G_{J}(b,0)=\frac{1}{2}F_{\theta_{J}}^{\prime\prime}(b)$

for any $b\in\mathbb{R}$, where $\partial_{2}G_{J}$ denotes the partial
derivative of $G_{J}$ with respect to its second argument. Computing this
derivative and using the two conditions yields

$\frac{1}{2}F_{\theta_{J}}^{\prime\prime}(b)=F_{\theta_{J}}^{\prime}(b)(\theta_{J}-b)\partial_{2}g_{J}(b,0),$

leading us to define

$g_{J}(b,v)\vcentcolon=\frac{v}{2(\theta_{J}-b)}\frac{F_{\theta_{J}}^{\prime\prime}(b)}{F_{\theta_{J}}^{\prime}(b)}$
(11)

for $b\neq\theta_{J}$, assuming that $F_{\theta_{J}}^{\prime}>0$. As long as
the distribution of the thresholds in population $J$ is unimodal and
symmetric, $F_{\theta_{J}}^{\prime\prime}(b)$ always has the same sign as
$\theta_{J}-b$, and in particular
$F_{\theta_{J}}^{\prime\prime}(\theta_{J})=0$. Hence, $g_{J}$ can be
continuously extended to $b=\theta_{J}$ by definition of the derivative
(provided $F_{\theta_{J}}^{\prime\prime\prime}(\theta_{J})$ exists), and it is
nonnegative on $\mathbb{R}\times[0,\infty)$.

For the approximation given by eq. (8) with $g_{J}$ defined as above to make
sense, we must ensure that $G_{J}$ has several properties. First,
$G_{J}\bigl{(}\E{B}[J](t),\C{BB}[JJ](t)\bigr{)}$ should still be consistent
with $\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]$ in extreme cases: if
$B_{t}^{J}$ has zero variance, then it should be possible to replace it with
its mean, while if it has infinite variance then $F_{\theta_{J}}(B_{t}^{J})$
should be $0$ or $1$ with $\nicefrac{{1}}{{2}}$–$\nicefrac{{1}}{{2}}$
probabilities. Additionally, if the expected input $\E{B}[J](t)$ goes to
$+\infty$ or $-\infty$, then $F_{\theta_{J}}(B_{t}^{J})$ should go to $1$ or
$0$ regardless of the variance. Moreover, since $F_{\theta_{J}}$ is
increasing, $G_{J}\bigl{(}\E{B}[J](t),\C{BB}[JJ](t)\bigr{)}$ should always
increase with $\E{B}[J](t)$ if the variance is fixed.

Adding some assumptions on the thresholds’ distribution allows to prove that
$G_{J}$ indeed has all of the properties enumerated above, in addition to
being consistent with the Taylor approximation of $F_{\theta_{J}}$ given in
eq. (10).

######  Theorem 1.

Suppose that the thresholds in population $J$ follow a unimodal and symmetric
distribution with mean $\theta_{J}$ and cumulative distribution function
$F_{\theta_{J}}$. Let $g_{J}$ and $G_{J}$ be defined by eq. (11) and eq. (9)
respectively. Suppose that

  1. 1.

$F_{\theta_{J}}$ is $\mathscr{C}^{4}$ on $\mathbb{R}$;

  2. 2.

$F_{\theta_{J}}^{\prime}>0$;

  3. 3.

$g_{J}(\cdot,1)$ is bounded on $\mathbb{R}$;

  4. 4.

$\forall b\in\mathbb{R}$, $g_{J}(b,1)+(\theta_{J}-b)\partial_{1}g_{J}(b,1)\geq
0$.

Then $G_{J}$ satisfies the following conditions.

  1. i.

$G_{J}$ is $\mathscr{C}^{1}$ on $\mathbb{R}\times[0,\infty)$;

  2. ii.

$G_{J}(\cdot,0)=F_{\theta_{J}}$ and $\forall
b\in\mathbb{R},G_{J}(b,v)\to\nicefrac{{1}}{{2}}$ as $v\to\infty$;

  3. iii.

$\forall v\geq 0,G_{J}(b,v)\to 0$ as $b\to-\infty$ and $G_{J}(b,v)\to 1$ as
$b\to\infty$;

  4. iv.

$\forall v\geq 0$, $G_{J}(\cdot,v)$ is increasing;

  5. v.

$\partial_{2}G_{J}(\cdot,0)=\frac{1}{2}F_{\theta_{J}}^{\prime\prime}$.

The proof is given in Appendix B. Assumption 4 does not have a simple
interpretation, but we use it to prove that the maps $G_{J}(\cdot,v)$ are
increasing, and it is fairly easy to verify for given distributions. In
particular, we show in Appendix B that if the thresholds follow either a
normal or a logistic distribution, then all assumptions hold. An example of
the function $G_{J}(\cdot,v)$ for a fixed variance $v$ is shown on Fig. 2.

### 3.2 The moment closure

Approximations to expectations of the form
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$ and
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]$ can be
constructed using the function $G_{J}$ from the previous section. We start by
giving an approximation of
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$ in terms of the
expectations $\E{S}[J](t)$ and $\E{B}[J](t)$ as well as the covariance
$\C{SB}[JJ](t)$ and the variance $\C{BB}[JJ](t)$. As in the simpler case of
the expectation of $F_{\theta_{J}}(B_{t}^{J})$, the approximation of
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$ must meet some
requirements. First, the difference between the approximation and
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})]\E{S}[J](t)$ should approximate
the covariance between $F_{\theta_{J}}(B_{t}^{J})$ and $S_{t}^{J}$, so its
sign should always follow that of the covariance $\C{SB}[JJ](t)$ because
$F_{\theta_{J}}$ is increasing. Then, if $\E{B}[J](t)$ goes to $+\infty$ or
$-\infty$, $F_{\theta_{J}}(B_{t}^{J})$ should go to either $1$ or $0$, so
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$ should go to either
$\E{S}[J](t)$ or $0$. Finally, by using a Taylor approximation of
$F_{\theta_{J}}$ and neglecting all third central moments, one can obtain the
approximation

$\displaystyle\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]$
$\displaystyle\approx F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}\E{S}[J](t)$
$\displaystyle\qquad+\frac{1}{2}F_{\theta_{J}}^{\prime\prime}\bigl{(}\E{B}[J](t)\bigr{)}\E{S}[J](t)\C{BB}[JJ](t)$
$\displaystyle\qquad+F_{\theta_{J}}^{\prime}\bigl{(}\E{B}[J](t)\bigr{)}\C{SB}[JJ](t),$

which should be valid when covariances are small. This gives conditions on the
value of the approximation and of its partial derivatives when
$\C{SB}[JJ](t)=\C{BB}[JJ](t)=0$.

Now, the function $G_{J}$ defined in eq. (9) can be used to construct the
desired approximation. The following result follows from Theorem 1 and is
proven in Appendix B.

######  Corollary 1.

Suppose that all assumptions of Theorem 1 hold. Then
$f\colon(0,\infty)\times\mathbb{R}^{3}\to(0,\infty)$ defined by

$f(s,b,c,v)\vcentcolon=sG_{J}\Bigl{(}b+\frac{c}{s},v\Bigr{)}$

satisfies the following conditions for any $s>0$, $b,c\in\mathbb{R}$ and
$v\geq 0$.

  1. i.

$f$ is $\mathscr{C}^{1}$ on $(0,\infty)\times\mathbb{R}^{3}$;

  2. ii.

$f(s,b,c,v)\gtreqless sG_{J}(b,v)$ when $c\gtreqless 0$, and in particular
$f(s,b,0,v)=sG_{J}(b,v)$ and $f(s,b,0,0)=sF_{\theta_{J}}(b)$;

  3. iii.

$f(s,b,c,v)\to 0$ as $b\to-\infty$ and $f(s,b,c,v)\to s$ as $b\to\infty$;

  4. iv.

$\partial_{3}f(s,b,0,0)=F_{\theta_{J}}^{\prime}(b)$;

  5. v.

$\partial_{4}f(s,b,0,0)=\frac{1}{2}sF_{\theta_{J}}^{\prime\prime}(b)$.

These results motivate the approximation

$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}]\approx\E{S}[J](t)G_{J}\Bigl{(}\E{B}[J](t)+\frac{\C{SB}[JJ](t)}{\E{S}[J](t)},\C{BB}[JJ](t)\Bigr{)}.$
(12)

Approximations to
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]$ and
$\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}R_{t}^{K}]$ can be
constructed in a similar way. Now, the only obvious conditions that the
approximation should satisfy are those related to the Taylor approximation of
$F_{\theta_{J}}$. Indeed, neglecting third central moments as before,

$\displaystyle\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]$
$\displaystyle\approx\bigl{(}\E{A}[K](t)\E{S}[J](t)+\C{AS}[KJ](t)\bigr{)}F_{\theta_{J}}\bigl{(}\E{B}[J](t)\bigr{)}$
$\displaystyle+\bigl{(}\E{A}[K](t)\C{SB}[JJ](t)+\E{S}[J](t)\C{AB}[KJ](t)\bigr{)}F_{\theta_{J}}^{\prime}\bigl{(}\E{B}[J](t)\bigr{)}$
$\displaystyle+\frac{1}{2}\E{A}[K](t)\E{S}[J](t)F_{\theta_{J}}^{\prime\prime}\bigl{(}\E{B}[J](t)\bigr{)}\C{BB}[JJ](t),$

which should hold when covariances are small. The appropriate approximation is
defined in the following result, which follows from Theorem 1 and is proven in
Appendix B.

######  Corollary 2.

Suppose that all assumptions of Theorem 1 hold. Then
$f\colon(0,\infty)^{2}\times\mathbb{R}^{5}\to\mathbb{R}$ defined by

$f(x,s,b,c_{1},c_{2},c_{3},v)\vcentcolon=(xs+c_{1})G_{J}\Bigl{(}b+\frac{c_{2}}{x}+\frac{c_{3}}{s},v\Bigr{)}$

satisfies the following conditions for any $x,s>0$, $b\in\mathbb{R}$ and
$v\geq 0$.

  1. i.

$f$ is $\mathscr{C}^{1}$ on $(0,\infty)^{2}\times\mathbb{R}^{5}$;

  2. ii.

$f(x,s,b,0,0,0,v)=xsG_{J}(b,v)$, and in particular
$f(x,s,b,0,0,0,0)=xsF_{\theta_{J}}(b)$;

  3. iii.

$\partial_{4}f(x,s,b,0,0,0,0)=F_{\theta_{J}}(b)$;

  4. iv.

$\partial_{5}f(x,s,b,0,0,0,0)=sF_{\theta_{J}}^{\prime}(b)$;

  5. v.

$\partial_{6}f(x,s,b,0,0,0,0)=xF_{\theta_{J}}^{\prime}(b)$;

  6. vi.

$\partial_{7}f(x,s,b,0,0,0,0)=\frac{1}{2}xsF_{\theta_{J}}^{\prime\prime}(b)$.

These results motivate the approximation

$\begin{multlined}\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]\approx\bigl{(}\E{A}[K](t)\E{S}[J](t)+\C{AS}[KJ](t)\bigr{)}\\\
\times
G_{J}\Bigl{(}\E{B}[J](t)+\frac{\C{AB}[KJ](t)}{\E{A}[K](t)}+\frac{\C{SB}[JJ](t)}{\E{S}[J](t)},\C{BB}[JJ](t)\Bigr{)}\end{multlined}\mathbb{E}^{\eta}[F_{\theta_{J}}(B_{t}^{J})S_{t}^{J}A_{t}^{K}]\approx\bigl{(}\E{A}[K](t)\E{S}[J](t)+\C{AS}[KJ](t)\bigr{)}\\\
\times
G_{J}\Bigl{(}\E{B}[J](t)+\frac{\C{AB}[KJ](t)}{\E{A}[K](t)}+\frac{\C{SB}[JJ](t)}{\E{S}[J](t)},\C{BB}[JJ](t)\Bigr{)}$
(13)

as well as a similar one where $A$ is replaced with $R$.

### 3.3 The dynamical system

The last approximation together with eq. (12) yield a closed dynamical system:

$\displaystyle\dE{A}[J]$
$\displaystyle=-\beta_{J}\E{A}[J]+\alpha_{J}\E{S}[J]G_{J}\Bigl{(}\E{B}[J]+\frac{\C{SB}[JJ]}{\E{S}[J]},\C{BB}[JJ]\Bigr{)},$
(14a) $\displaystyle\dE{R}[J]$
$\displaystyle=-\gamma_{J}\E{R}[J]+\beta_{J}\E{A}[J],$ (14b)
$\displaystyle\begin{split}\dC{AA}[JK]&=-(\beta_{J}+\beta_{K})\C{AA}[JK]\\\
&\quad+\alpha_{K}H_{K}(\E{A}[J],\E{S}[K],\E{B}[K],\C{AS}[JK],\C{AB}[JK],\C{SB}[KK],\C{BB}[KK])\\\
&\quad+\alpha_{J}H_{J}(\E{A}[K],\E{S}[J],\E{B}[J],\C{AS}[KJ],\C{AB}[KJ],\C{SB}[JJ],\C{BB}[JJ]),\end{split}$
(14c) $\displaystyle\dC{RR}[JK]$
$\displaystyle=-(\gamma_{J}+\gamma_{K})\C{RR}[JK]+\beta_{K}\C{AR}[KJ]+\beta_{J}\C{AR}[JK],$
(14d)
$\displaystyle\begin{split}\dC{AR}[JK]&=-(\beta_{J}+\gamma_{K})\C{AR}[JK]+\beta_{K}\C{AA}[JK]\\\
&\quad+\alpha_{J}H_{J}(\E{R}[K],\E{S}[J],\E{B}[J],\C{RS}[KJ],\C{RB}[KJ],\C{SB}[JJ],\C{BB}[JJ]),\end{split}$
(14e)

where

$H_{J}(x,s,b,c_{1},c_{2},c_{3},v)\vcentcolon=(xs+c_{1})G_{J}\Bigl{(}b+\frac{c_{2}}{x}+\frac{c_{3}}{s},v\Bigr{)}\\\
-xsG_{J}\Bigl{(}b+\frac{c_{3}}{s},v\Bigr{)}.$

For a network of $n$ populations, there are $2n$ equations for expectations of
active and refractory fractions, $n(n+1)$ equations for covariances between
two active or two refractory fractions (including variances), and $n^{2}$
equations for covariances between an active and a refractory fraction. Hence,
the whole system has $n(2n+3)$ dimensions.

This dynamical system can be interpreted from the definition of the function
$G_{J}$. In particular, eq. (14a) has the same form as the corresponding
equation in the mean-field system, but the activation function is stretched
horizontally by the variance of the input and shifted by the covariance
between the input and the sensitive fraction. The derivatives of covariances
have similar interpretations.

Corollaries to Theorem 1 guarantee that the functions $G_{J}$ and $H_{J}$ are
always $\mathscr{C}^{1}$ on the relevant domains. Therefore, the vector field
corresponding to eq. (14) is $\mathscr{C}^{1}$ on
$(0,1)^{2n}\times\mathbb{R}^{n(2n+1)}$. The Picard–Lindelöf theorem then shows
the existence and uniqueness of solutions to the differential equation for any
initial condition in $(0,1)^{2n}\times\mathbb{R}^{n(2n+1)}$. However, not
every initial condition in this domain can be interpreted in terms of the
underlying Markov chain, since the boundedness of the fractions of populations
implies various bounds on the covariances, for example through the
Bhatia–Davis inequality [4] or the Cauchy–Schwarz inequality. Given the
various conditions that relate all covariances together, it is not obvious to
usefully characterize the domain where solutions certainly make sense from the
point of view of the underlying Markov chain, but we do not expect this to be
a problem as long as covariances remain small.

It is worth noting the following result about the case where all covariances
are zero.

######  Proposition 1.

The domain $\mathscr{D}\vcentcolon=(0,1)^{2n}\times\\{0\\}^{n(2n+1)}$, where
$\C{AA}[JK]=\C{RR}[JK]=\C{AR}[JK]=0$ for all $J,K\in\mathscr{P}$, is invariant
under the flow of the dynamical system given in eq. (14).

###### Proof.

If $\C{AA}[JK]=\C{RR}[JK]=\C{AR}[JK]=0$ for all $J,K\in\mathscr{P}$, then all
other covariances are zero as well. It is then clear from the definition of
the functions $H_{J}$ and from eq. (14) that
$\dC{AA}[JK]=\dC{RR}[JK]=\dC{AR}[JK]=0$. Hence, on the domain $\mathscr{D}$,
the vector field corresponding to the dynamical system is parallel to
$\mathscr{D}$, which implies its invariance. ∎

Proposition 1 has an important consequence, since by Theorem 1,
$G_{J}(\E{B}[J],0)=F_{\theta_{J}}(\E{B}[J])$. Indeed, it follows that if all
covariances are forced to zero, the system reduces to the mean-field system
given in eq. (6). This shows that we can see the mean-field system as a
subsystem of that given in eq. (14). Since this mean-field system can be seen
as an extension of Wilson–Cowan’s system [48], this means that our second-
order system is also an extension of Wilson–Cowan’s.

## 4 Examples

We now illustrate how adding covariances can lead to a better prediction of
the underlying Markov chain’s macroscopic behavior. We first give an example
where the second-order model correctly predicts the steady state to which
trajectories of the underlying Markov chain converge, whereas the mean-field
model does not. Then, we present two examples in which the second-order model
converges to a steady state with nonzero covariances that cannot be described
by the mean-field model. In both cases, while the mean-field model is able to
predict some aspects of the Markov chain’s macroscopic dynamics, the second-
order model provides more information, leading to a better prediction of the
Markov chain’s behavior.

Figure 3: Numerical simulations for a network of a single population of 1000
neurons with parameters given in eq. (15) and initial state given in eq. (16).
On the left panels is a comparison between the solution of the second-order
dynamical system and statistics computed from 1000 simulated trajectories of
the underlying Markov chain. The shaded regions around the curves are
associated with statistics, and are bounded above and below by a difference of
one standard deviation from the average value. On the right panel is a phase
plane of the mean-field dynamical system, on which are plotted the solution of
the mean-field system and the macroscopic behavior of a typical trajectory of
the Markov chain.

### 4.1 A word on the methodology

In all examples, we compare the solutions of the mean-field and second-order
models with trajectories of the underlying Markov chain obtained using the
Doob–Gillespie algorithm [27, 28]. To produce a meaningful comparison between
the results of macroscopic and microscopic models, it is crucial to ensure
that the microscopic parameters and initial states are consistent with the
macroscopic ones.

For simplicity, we take the microscopic parameters $\alpha$, $\beta$ and
$\gamma$ constant over populations and equal to their macroscopic values. In
the same way, having chosen a connection coefficient $c_{JK}$ from a
population $K$ to a population $J$, we set the weight of the connection from
any neuron $k\in K$ to any $j\in J$ to be $W_{jk}=\nicefrac{{c_{JK}}}{{\lvert
K\rvert}}$. However, we do not assume the thresholds to be constant: we rather
assume that thresholds in a population $J$ follow a logistic distribution with
mean $\theta_{J}$ and scaling factor $s_{\theta_{J}}$, having cumulative
distribution function

$F_{\theta_{J}}(x)=\frac{1}{1+\exp\Bigl{(}-\dfrac{x-\theta_{J}}{s_{\theta_{J}}}\Bigr{)}}.$

Finally, we always specify the distribution of the microscopic initial state
from given values for the initial expectations of the active and refractory
fractions of each population. To do this, we assume that the initial states of
all neurons are independent and that they are identically distributed over
populations. Indeed, given initial expectations $\E{A}[J](0),\E{R}[J](0)$ for
a population $J$, this assumption implies that the state of each neuron $j\in
J$ is $1$, $i$ or $0$ with probabilities $\E{A}[J](0)$, $\E{R}[J](0)$ or
$\E{S}[J](0)$ respectively. Then, the values of all initial covariances can
easily be related to the initial expectations. For example,

$\C{AA}[JJ](0)=\frac{1}{\lvert J\rvert^{2}}\sum_{j\in
J}\mathrm{Var}\big{[}\operatorname{Re}X_{j}^{0}\big{]}=\frac{1}{\lvert
J\rvert}\E{A}[J](0)\bigl{(}1-\E{A}[J](0)\bigr{)}.$

The same pattern works for the other neuronal states, and $\C{AR}[JJ](0)$ can
be computed from the initial variances. Since the initial states of neurons
are independent, the initial covariance between fractions of two distinct
populations is zero. Using this method, all initial covariances are determined
by the initial expectations once the size of each population is fixed.

### 4.2 Covariances can improve the model’s accuracy

The first example is a case where the second-order model allows to predict
accurately the macroscopic behavior of the underlying Markov chain, while the
mean-field model does not.

Consider a network of $N$ neurons in a single population, with parameters

$\displaystyle\alpha$ $\displaystyle=1.4\,[\gamma],$ $\displaystyle\beta$
$\displaystyle=2.5\,[\gamma],$ $\displaystyle\gamma$
$\displaystyle=1\,[\gamma],$ (15) $\displaystyle\theta$ $\displaystyle=0.75,$
$\displaystyle s_{\theta}$ $\displaystyle=0.1,$ $\displaystyle Q$
$\displaystyle=0,$ $\displaystyle c$ $\displaystyle=5.5,$

where we dropped the subscripts that would refer to the unique population.
Here, we measure characteristic rates in units of $\gamma$, which is
equivalent to measuring time in units of $\nicefrac{{1}}{{\gamma}}$ since
every term of the dynamical system is proportional to one of the rates. We fix
initial expected values

$(\E{A},\E{R})(0)=(0.16,0.51),$ (16)

which determines the distribution of the microscopic initial state (and thus
the macroscopic initial covariances) as a function of $N$. Here, we use a
network of $N=1000$ neurons.

Integrating numerically the second-order system with parameters from eq. (15)
and initial state given by eq. (16) yields the solution shown on the left
panels of Fig. 3. However, integrating the mean-field system with the same
parameters and initial state yields the solution shown on the phase plane on
the right panel of Fig. 3. The two models disagree: the mean-field model
predicts that the system converges to a fixed point with about $20\%$ of
active neurons, while the second-order model predicts that all activity will
stop.

It is interesting to compare these solutions to statistics computed from
trajectories of the underlying Markov chain. Choosing randomly microscopic
initial states corresponding to the macroscopic initial state given by eq.
(16), we computed the averages and sample covariances of the three fractions
of the network in each of the three states over 1000 simulated trajectories.
The results are shown on the left panels of Fig. 3. In this case, the second-
order model predicts more accurately the macroscopic behavior of the
underlying Markov chain than the mean-field model.

### 4.3 The second-order model can describe the averaging of multiple steady
states

Figure 4: Numerical simulations for a network consisting of a single
population of 100 neurons with parameters given in eq. (17) and initial state
given in eq. (18). On the left panels is a comparison between the solution of
the second-order dynamical system and statistics computed from 10 000
simulated trajectories of the underlying Markov chain. On the right panel is a
phase plane of the mean-field dynamical system, on which are plotted the
solution of the mean-field system and the macroscopic behaviors of two
trajectories of the Markov chain. These trajectories are typical
representatives of the sets of trajectories that converge to each fixed point.

This example illustrates an aspect of the macroscopic dynamics of the Markov
chain that can be reproduced by our second-order model and not by the mean-
field model.

Consider a network of $N$ neurons in a single population, with parameters

$\displaystyle\alpha$ $\displaystyle=4.2\,[\gamma],$ $\displaystyle\beta$
$\displaystyle=0.05\,[\gamma],$ $\displaystyle\gamma$
$\displaystyle=1\,[\gamma],$ (17) $\displaystyle\theta$ $\displaystyle=12.7,$
$\displaystyle s_{\theta}$ $\displaystyle=0.2,$ $\displaystyle Q$
$\displaystyle=0,$ $\displaystyle c$ $\displaystyle=17.$

The phase plane of the mean-field system on the right panel of Fig. 4 shows
that this system has two stable fixed points, one inactive at the origin, and
the other very active, with about 95% of active neurons. Now, fix the initial
expectations

$(\E{A},\E{R})(0)=(0.71,0.221).$ (18)

From this state, integrating the mean-field system yields the solution shown
on the phase plane, which converges to the inactive fixed point.

The initial condition given by eq. (18) is very close to the separatrix
between the basins of attraction of the two fixed points. As one can expect,
it is verified by stochastic simulations of the Markov chain that when
starting from a microscopic initial state drawn from the distribution
specified by eq. (18), noise can make trajectories converge to either of the
two stable fixed points. Two examples of trajectories are shown on the phase
plane on Fig. 4. The second-order model can predict this behavior.

Integrating the second-order system with initial state given by eq. (18) for a
network of $N=100$ neurons yields the solution shown on the left panels of
Fig. 4. It can be seen that the state to which the solution converges is the
average of the two stable fixed points of the mean-field system. Moreover, the
variances of the active and sensitive fractions of the network are close to
$\nicefrac{{1}}{{4}}$, which is seen from the Bhatia–Davis inequality [4] to
be the upper bound on the variance of random variables supported on $[0,1]$.
To interpret what it means for the variances of $A_{t}$ and $S_{t}$ to have
such high values, recall that the maximal variance of $\nicefrac{{1}}{{4}}$ is
attained by a random variable that is $0$ or $1$ with probabilities
$\nicefrac{{1}}{{2}}$. Thus, the results shown on the left panels of Fig. 4
can represent that when $t$ is large enough, $A_{t}$ and $S_{t}$ are close to
binary random variables with values at the two fixed points of the mean-field
system.

To compare this solution with the macroscopic behavior of the underlying
Markov chain, we generated 10 000 trajectories of the stochastic process with
a network of $N=100$ neurons. The relevant statistics from these trajectories
are shown on the left panels of Fig. 4. The results of the stochastic
simulations agree well with the prediction of the second-order model. On the
other hand, the mean-field model is only a good approximation of the
trajectories that converge to the inactive fixed point—it fails to capture the
other possible outcome.

### 4.4 The second-order model can describe the averaging of oscillations

This last example illustrates again the ability of the second-order model to
carry information about the distribution of the underlying Markov chain that
is inaccessible to the mean-field model.

Consider a network split into two populations $E$ and $I$, the former being
excitatory with parameters

$\displaystyle\alpha_{E}$ $\displaystyle=0.75\,[\gamma_{E}],$
$\displaystyle\beta_{E}$ $\displaystyle=0.15\,[\gamma_{E}],$
$\displaystyle\gamma_{E}$ $\displaystyle=1\,[\gamma_{E}],$ (19a)
$\displaystyle\theta_{E}$ $\displaystyle=0.7,$ $\displaystyle s_{\theta_{E}}$
$\displaystyle=0.2,$ $\displaystyle Q_{E}$ $\displaystyle=0,$ (19b) and the
latter being inhibitory with parameters $\displaystyle\alpha_{I}$
$\displaystyle=0.4\,[\gamma_{E}],$ $\displaystyle\beta_{I}$
$\displaystyle=0.12\,[\gamma_{E}],$ $\displaystyle\gamma_{I}$
$\displaystyle=0.5\,[\gamma_{E}],$ (19c) $\displaystyle\theta_{I}$
$\displaystyle=1.8,$ $\displaystyle s_{\theta_{I}}$ $\displaystyle=0.2,$
$\displaystyle Q_{I}$ $\displaystyle=0.$ (19d) The connections between these
populations are given by $c=\begin{pmatrix}c_{EE}&c_{EI}\\\
c_{IE}&c_{II}\end{pmatrix}=\begin{pmatrix}11&-12\\\ 12&-9\end{pmatrix}.$ (19e)

Now, fix the initial state from the expectations

$(\E{A}[E],\E{A}[I],\E{R}[E],\E{R}[I])(0)=(0.25,0.3,0.2,0.25).$ (20)

The solution of the mean-field system converges to a limit cycle, as shown on
the top left panel of Fig. 5. This solution approximates reasonably well the
behavior of the underlying Markov chain, as seen from the example of
trajectory shown on the top right panel of Fig. 5.

Figure 5: Numerical simulations for a network of two populations with
parameters given in eq. (19) and initial state given in eq. (20), with 100
neurons in each population. On the top left panel is the solution of the mean-
field dynamical system. On the top right panel is the macroscopic behavior of
a typical trajectory of the underlying Markov chain obtained from a numerical
simulation. On the middle and bottom panels is a comparison between the
solution of the second-order system (left panels) and statistics computed from
10 000 simulated trajectories of the Markov chain (right panels).

Yet, one can wonder what happens when statistics are computed from the
macroscopic behavior of multiple trajectories. As the trajectories are random,
the oscillation periods are expected to vary among trajectories, and even
among oscillations in each trajectory. Thus, one should expect the
oscillations to decay over time when averaging many trajectories. This is the
behavior that our second-order model reproduces.

Integrating the second-order system with initial state given by eq. (20) for
populations of 100 neurons each yields the solution shown on the middle and
bottom left panels of Fig. 5. The solution converges to a fixed point where
variances are nonzero, and where the value of each expectation is about the
average, over a cycle, of the corresponding expectation predicted by the mean-
field system. To compare this solution with the macroscopic behavior of the
underlying Markov chain, we generated 10 000 trajectories in a network with
100 neurons in each population. The relevant statistics from these
trajectories are shown on the middle and bottom right panels of Fig. 5 and
illustrate that, at least qualitatively, the second-order model succeeds in
predicting the macroscopic behavior of the Markov chain.

## 5 Discussion

We demonstrated the importance of accounting for covariances between states of
neural populations in large stochastic neural networks. Our observations were
based on the low-dimensional dynamical system in eq. (14), which we derived
mostly by developing a novel method for estimating the expected value of
random variables as in eq. (12) and eq. (13), and solving the corresponding
moment closure problem. Given the nonlinearity of the model and the presence
of multiple attractors, obtaining rigorous theorems on rates of convergence or
on the quality of our dimension reduction would be difficult. Nevertheless,
our numerical experiments showed that in some cases, the mean-field solutions
greatly differ from both the average of realizations of the stochastic model
and the solutions of the second-order model.

There is a first class of examples (see Fig. 3 and Fig. 4) in which the
inclusion of covariances leads to predictions that match much more closely the
behavior of the high-dimensional stochastic system. This can occur even if
after a brief transient period, all covariances go to zero. The presence of
two attractors in the mean-field system explains well this phenomenon. Indeed,
if the initial conditions lie close to the separatrix between the basins of
attraction, the qualitative features of the solution can be highly sensitive
to small perturbations of the early system behavior. Failing to account for
covariances is akin to not using all the available information about the
system’s initial state, which may lead to predicting convergence to the wrong
attractor. A lot of research has been done recently for systems with only one
attracting point [25, 24, 23]. In this case, injecting information about
covariances increases the order of convergence. In our case, covariances have
more dramatic impacts as they may make the system jump basin of attraction.

There is yet another situation in which the covariances may dramatically
affect the solutions as some covariances can remain positive throughout the
simulations. As shown on Fig. 5, the solution of the mean-field system can
display sustained oscillations while the solution of the second-order model
converges to a fixed point. This can be explained by the presence of positive
covariances, which impacts the activation functions of neuron populations by
making them more linear. Since the presence of multiple fixed points or
oscillations in the mean-field model depends on the nonlinearity of the
activation function, changing its shape is likely to change the topology of
the solutions.

We argue for the necessity of taking into account the inherently stochastic
nature of neuron activity. This stochastic nature does not come from the way
that mathematical modeling of neurons is performed, but from omnipresent
biological mechanisms such as vesicle release or random channel state
transition. Our results show that using a naive approach to handle this
stochasticity, such as a mean-field approximation, can lead to false
conclusions. Our work thus emphasizes the importance of pursuing the study of
higher moments in models of neuronal activity. In particular, the observation
that the solutions of the second-order model correspond to the averages of
many trajectories of the stochastic model invites further theoretical
investigation. We believe that the mathematical tools developed here will help
deriving more effective low-dimensional approximations of stochastic processes
on networks, contributing to a better understanding of large-scale dynamical
phenomena in complex systems.

#### Code Availability

All numerical results presented in this paper were obtained with the PopNet
package [47], written in the Python programming language and available on
GitHub.

#### Acknowledgments

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (N.D, P.D, V.P), the Fonds de recherche du Québec – Nature
and technologies (N.D., V.P.), and the Sentinel North program of Université
Laval (N.D., P.D), funded by the Canada First Research Excellence Fund.

## References

  * [1] Bruno B. Averbeck, Peter E. Latham and Alexandre Pouget  “Neural Correlations, Population Coding and Computation”  In _Nature Reviews Neuroscience_ 7.5 Nature Publishing Group, 2006, pp. 358–366  DOI: 10.1038/nrn1888
  * [2] Rava Azeredo da Silveira and Fred Rieke  “The Geometry of Information Coding in Correlated Neural Populations”  In _Annual Review of Neuroscience_ 44.1, 2021, pp. 403–424  DOI: 10.1146/annurev-neuro-120320-082744
  * [3] R.. Beurle  “Properties of a Mass of Cells Capable of Regenerating Pulses”  In _Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences_ 240.669 Royal Society, 1956, pp. 55–94  DOI: 10.2307/92540
  * [4] Rajendra Bhatia and Chandler Davis  “A Better Bound on the Variance”  In _The American Mathematical Monthly_ 107.4 Taylor & Francis, 2000, pp. 353–357  DOI: 10.1080/00029890.2000.12005203
  * [5] Paul C. Bressloff  “Path-Integral Methods for Analyzing the Effects of Fluctuations in Stochastic Hybrid Neural Networks”  In _The Journal of Mathematical Neuroscience (JMN)_ 5.1, 2015, pp. 4  DOI: 10.1186/s13408-014-0016-z
  * [6] Paul C. Bressloff  “Stochastic Neural Field Theory and the System-Size Expansion”  In _SIAM Journal on Applied Mathematics_ 70.5 Society for Industrial and Applied Mathematics, 2009, pp. 1488–1521  DOI: 10.1137/090756971
  * [7] Paul C. Bressloff, G. Ermentrout, Olivier Faugeras and Peter J. Thomas  “Stochastic Network Models in Neuroscience: A Festschrift for Jack Cowan. Introduction to the Special Issue”  In _Journal of Mathematical Neuroscience_ 6.1 Springer Verlag, 2016, pp. 4  DOI: 10.1186/s13408-016-0036-y
  * [8] Nicolas Brunel  “Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons”  In _Journal of Computational Neuroscience_ 8.3, 2000, pp. 183–208  DOI: 10.1023/A:1008925309027
  * [9] Michael A. Buice and Carson C. Chow  “Beyond Mean Field Theory: Statistical Field Theory for Neural Networks”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2013.03 IOP Publishing, 2013, pp. P03003  DOI: 10.1088/1742-5468/2013/03/P03003
  * [10] Michael A. Buice and Jack D. Cowan  “Field-Theoretic Approach to Fluctuation Effects in Neural Networks”  In _Physical Review E_ 75.5 American Physical Society, 2007, pp. 051919  DOI: 10.1103/PhysRevE.75.051919
  * [11] Michael A. Buice and Jack D. Cowan  “Statistical Mechanics of the Neocortex”  In _Progress in Biophysics and Molecular Biology_ 99.2, 2009, pp. 53–86  DOI: 10.1016/j.pbiomolbio.2009.07.003
  * [12] Michael A. Buice, Jack D. Cowan and Carson C. Chow  “Systematic Fluctuation Expansion for Neural Network Activity Equations”  In _Neural Computation_ 22.2 MIT Press, 2010, pp. 377–426  DOI: 10.1162/neco.2009.02-09-960
  * [13] W H Calvin and C F Stevens  “Synaptic Noise and Other Sources of Randomness in Motoneuron Interspike Intervals.”  In _Journal of Neurophysiology_ 31.4 American Physiological Society, 1968, pp. 574–587  DOI: 10.1152/jn.1968.31.4.574
  * [14] E. Cator and P. Van Mieghem  “Second-Order Mean-Field Susceptible-Infected-Susceptible Epidemic Threshold”  In _Physical Review E_ 85.5 American Physical Society, 2012, pp. 056111  DOI: 10.1103/PhysRevE.85.056111
  * [15] Carson C. Chow and Yahya Karimipanah  “Before and beyond the Wilson–Cowan Equations”  In _Journal of Neurophysiology_ 123.5 American Physiological Society, 2020, pp. 1645–1656  DOI: 10.1152/jn.00404.2019
  * [16] Jack D. Cowan  “Stochastic Neurodynamics”  In _Advances in Neural Information Processing Systems_ 3 Morgan-Kaufmann, 1990, pp. 62–69  URL: https://proceedings.neurips.cc/paper/1990/hash/3c7781a36bcd6cf08c11a970fbe0e2a6-Abstract.html
  * [17] Jack D. Cowan, Jeremy Neuman and Wim Drongelen  “Wilson–Cowan Equations for Neocortical Dynamics”  In _The Journal of Mathematical Neuroscience_ 6.1 SpringerOpen, 2016, pp. 1  DOI: 10.1186/s13408-015-0034-5
  * [18] Alain Destexhe and Terrence J Sejnowski  “The Wilson-Cowan Model, 36 Years Later”  In _Biological Cybernetics_ 101, 2009, pp. 1–2  DOI: 10.1007/s00422-009-0328-3
  * [19] C. Domb  “On the Theory of Cooperative Phenomena in Crystals”  In _Advances in Physics_ 9.35 Taylor & Francis, 1960, pp. 245–361  DOI: 10.1080/00018736000101199
  * [20] Joseph L. Doob  “Stochastic Processes”, Wiley Classics Library  John Wiley & Sons, 1990
  * [21] G. Ermentrout and David H. Terman  “Mathematical Foundations of Neuroscience” 35, Interdisciplinary Applied Mathematics  New York, NY: Springer, 2010  DOI: 10.1007/978-0-387-87708-2
  * [22] A. Faisal, Luc P.. Selen and Daniel M. Wolpert  “Noise in the Nervous System”  In _Nature Reviews Neuroscience_ 9.4 Nature Publishing Group, 2008, pp. 292–303  DOI: 10.1038/nrn2258
  * [23] Nicolas Gast  “Refinements of Mean Field Approximation”, 2020  URL: https://theses.hal.science/tel-02509756
  * [24] Nicolas Gast, Luca Bortolussi and Mirco Tribastone  “Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis”  In _Performance Evaluation_ 129, 2019, pp. 60–80  DOI: 10.1016/j.peva.2018.09.005
  * [25] Nicolas Gast and Benny Van Houdt  “A Refined Mean Field Approximation”  In _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ 1.2, 2017, pp. 33  DOI: 10.1145/3154491
  * [26] Wulfram Gerstner  “Time Structure of the Activity in Neural Network Models”  In _Physical Review E_ 51.1 American Physical Society, 1995, pp. 738–758  DOI: 10.1103/PhysRevE.51.738
  * [27] Daniel T. Gillespie  “A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions”  In _Journal of Computational Physics_ 22, 1976, pp. 403–434  DOI: 10.1016/0021-9991(76)90041-3
  * [28] Daniel T. Gillespie  “Exact Stochastic Simulation of Coupled Chemical Reactions”  In _Journal of Physical Chemistry_ 81.25 American Chemical Society, 1977, pp. 2340–2361  DOI: 10.1021/j100540a008
  * [29] J.. Griffith  “A Field Theory of Neural Nets: I: Derivation of Field Equations”  In _The Bulletin of Mathematical Biophysics_ 25.1 Springer, 1963, pp. 111–120  DOI: 10.1007/BF02477774
  * [30] J.. Griffith  “A Field Theory of Neural Nets: II. Properties of the Field Equations”  In _The Bulletin of Mathematical Biophysics_ 27.2 Springer, 1965, pp. 187–195  DOI: 10.1007/BF02498774
  * [31] Haiping Huang  “Statistical Mechanics of Neural Networks”  Singapore: Springer Nature, 2021  DOI: 10.1007/978-981-16-7570-6
  * [32] Jaewook Joo and Joel L. Lebowitz  “Pair Approximation of the Stochastic Susceptible-Infected-Recovered-Susceptible Epidemic Model on the Hypercubic Lattice”  In _Physical Review E_ 70.3 American Physical Society, 2004, pp. 036114  DOI: 10.1103/PhysRevE.70.036114
  * [33] Stephen Keeley, Áine Byrne, André Fenton and John Rinzel  “Firing Rate Models for Gamma Oscillations”  In _Journal of Neurophysiology_ 121.6 American Physiological Society, 2019, pp. 2181–2190  DOI: 10.1152/jn.00741.2018
  * [34] M.. Keeling  “The Effects of Local Spatial Structure on Epidemiological Invasions”  In _Proceedings of the Royal Society of London. Series B: Biological Sciences_ 266.1421 Royal Society, 1999, pp. 859–867  DOI: 10.1098/rspb.1999.0716
  * [35] Istvan Z. Kiss et al.  “Exact Deterministic Representation of Markovian SIR Epidemics on Networks with and without Loops”  In _Journal of Mathematical Biology_ 70.3, 2015, pp. 437–464  DOI: 10.1007/s00285-014-0772-0
  * [36] Christian Kuehn  “Moment Closure—A Brief Review”  In _Control of Self-Organizing Nonlinear Systems_ , Understanding Complex Systems  Cham: Springer International Publishing, 2016, pp. 253–271  DOI: 10.1007/978-3-319-28028-8˙13
  * [37] Lior Lebovich et al.  “Idiosyncratic Choice Bias Naturally Emerges from Intrinsic Stochasticity in Neuronal Dynamics”  In _Nature Human Behaviour_ 3.11 Nature Publishing Group, 2019, pp. 1190–1202  DOI: 10.1038/s41562-019-0682-7
  * [38] J L LeBowitz  “Statistical Mechanics-A Review of Selected Rigorous Results”  In _Annual Review of Physical Chemistry_ 19.1, 1968, pp. 389–418  DOI: 10.1146/annurev.pc.19.100168.002133
  * [39] Simon A. Levin and Richard Durrett  “From Individuals to Epidemics”  In _Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences_ 351.1347 Royal Society, 1996, pp. 1615–1621  DOI: 10.1098/rstb.1996.0145
  * [40] David T.. Liley, Peter J. Cadusch and Mathew P. Dafilis  “A Spatially Continuous Mean Field Theory of Electrocortical Activity”  In _Network: Computation in Neural Systems_ 13.1, 2001, pp. 67  DOI: 10.1088/0954-898X/13/1/303
  * [41] Wolfgang Maass  “Noise as a Resource for Computation and Learning in Networks of Spiking Neurons”  In _Proceedings of the IEEE_ 102.5, 2014, pp. 860–880  DOI: 10.1109/JPROC.2014.2310593
  * [42] Hirotsugu Matsuda, Naofumi Ogita, Akira Sasaki and Kazunori Satō  “Statistical Mechanics of Population: The Lattice Lotka-Volterra Model”  In _Progress of Theoretical Physics_ 88.6, 1992, pp. 1035–1049  DOI: 10.1143/ptp/88.6.1035
  * [43] Max Nolte et al.  “Cortical Reliability amid Noise and Chaos”  In _Nature Communications_ 10.1 Nature Publishing Group, 2019, pp. 3792  DOI: 10.1038/s41467-019-11633-8
  * [44] James R. Norris  “Markov Chains”, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge: Cambridge University Press, 1997  DOI: 10.1017/CBO9780511810633
  * [45] Gabriel Koch Ocker, Krešimir Josić, Eric Shea-Brown and Michael A. Buice  “Linking Structure and Activity in Nonlinear Spiking Networks”  In _PLOS Computational Biology_ 13.6 Public Library of Science, 2017, pp. e1005583  DOI: 10.1371/journal.pcbi.1005583
  * [46] Vincent Painchaud  “Dynamique markovienne ternaire cyclique sur graphes et quelques applications en biologie mathématique”, 2021  URL: https://corpus.ulaval.ca/jspui/handle/20.500.11794/71616
  * [47] Vincent Painchaud  “PopNet”, 2022  DOI: 10.5281/zenodo.7458680
  * [48] Vincent Painchaud, Nicolas Doyon and Patrick Desrosiers  “Beyond Wilson–Cowan Dynamics: Oscillations and Chaos without Inhibition”  In _Biological Cybernetics_ 116.5, 2022, pp. 527–543  DOI: 10.1007/s00422-022-00941-w
  * [49] Stefano Panzeri, Monica Moroni, Houman Safaai and Christopher D. Harvey  “The Structures and Functions of Correlations in Neural Population Codes”  In _Nature Reviews Neuroscience_ 23.9 Nature Publishing Group, 2022, pp. 551–567  DOI: 10.1038/s41583-022-00606-4
  * [50] Emilio Salinas and Terrence J. Sejnowski  “Correlated Neuronal Activity and the Flow of Neural Information”  In _Nature Reviews Neuroscience_ 2.8 Nature Publishing Group, 2001, pp. 539–550  DOI: 10.1038/35086012
  * [51] Kazunori Satō, Hirotsugu Matsuda and Akira Sasaki  “Pathogen Invasion and Host Extinction in Lattice Structured Populations”  In _Journal of Mathematical Biology_ 32.3, 1994, pp. 251–268  DOI: 10.1007/BF00163881
  * [52] Sebastian Schneegans, Robert Taylor and Paul M. Bays  “Stochastic Sampling Provides a Unifying Account of Visual Working Memory Limits”  In _Proceedings of the National Academy of Sciences_ 117.34 Proceedings of the National Academy of Sciences, 2020, pp. 20959–20968  DOI: 10.1073/pnas.2004306117
  * [53] Elad Schneidman, Michael J. Berry, Ronen Segev and William Bialek  “Weak Pairwise Correlations Imply Strongly Correlated Network States in a Neural Population”  In _Nature_ 440.7087 Nature Publishing Group, 2006, pp. 1007–1012  DOI: 10.1038/nature04701
  * [54] Kieran J. Sharkey  “Deterministic Epidemiological Models at the Individual Level”  In _Journal of Mathematical Biology_ 57.3, 2008, pp. 311–331  DOI: 10.1007/s00285-008-0161-7
  * [55] Tim P. Vogels, Kanaka Rajan and L.F. Abbott  “Neural Network Dynamics”  In _Annual Review of Neuroscience_ 28.1, 2005, pp. 357–376  DOI: 10.1146/annurev.neuro.28.061604.135637
  * [56] John A. White, Jay T. Rubinstein and Alan R. Kay  “Channel Noise in Neurons”  In _Trends in Neurosciences_ 23.3, 2000, pp. 131–137  DOI: 10.1016/S0166-2236(99)01521-0
  * [57] Hugh R. Wilson and Jack D. Cowan  “Evolution of the Wilson–Cowan Equations”  In _Biological Cybernetics_ 115.6, 2021, pp. 643–653  DOI: 10.1007/s00422-021-00912-7
  * [58] Hugh R. Wilson and Jack D. Cowan  “Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons”  In _Biophysical Journal_ 12.1 Cell Press, 1972, pp. 1–24  DOI: 10.1016/S0006-3495(72)86068-5
  * [59] Mahdi Zarepour et al.  “Universal and Nonuniversal Neural Dynamics on Small World Connectomes: A Finite-Size Scaling Analysis”  In _Physical Review E_ 100.5 American Physical Society, 2019, pp. 52138  DOI: 10.1103/PhysRevE.100.052138
  * [60] Yann Zerlaut and Alain Destexhe  “Enhanced Responsiveness and Low-Level Awareness in Stochastic Network States”  In _Neuron_ 94.5, 2017, pp. 1002–1009  DOI: 10.1016/j.neuron.2017.04.001

##  Appendix A Derivation of the nonautonomous microscopic and macroscopic
evolution equations

In this section, we derive the nonautonomous evolution equations given in the
main text in eq. (4), eq. (5) and eq. (7). In order to derive eq. (7), we
obtain as an intermediate step yet another system that can be seen as a
second-order version of that given in eq. (4). We derive each of these systems
in the first four subsections, and then we conclude the section by some
remarks on the limits of the approximations that we use to obtain the
macroscopic systems.

### A.1 The nonautonomous microscopic first-order evolution equation

First, we derive the expressions of the derivatives of the probabilities

$p_{j}(t)\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=1]=\mathbb{E}^{\eta}\big{[}\operatorname{Re}X_{t}^{j}\big{]},\qquad
r_{j}(t)\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=i]=\mathbb{E}^{\eta}\big{[}\operatorname{Im}X_{t}^{j}\big{]}\qquad\text{and}\qquad
q_{j}(t)\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=0]=\mathbb{E}^{\eta}\big{[}1-\lvert
X_{t}^{j}\rvert\big{]}$

where, as in the main text, the continuous-time Markov chain
$\\{X_{t}\\}_{t\geq 0}$ takes its values in $E\vcentcolon=\\{0,1,i\\}^{N}$,
$N$ being the number of neurons in the network. An entry $X_{t}^{j}$ is then
interpreted as the state of neuron $j$ at time $t$.

To obtain the expression of the derivative of $p_{j}(t)$, let $\Delta t>0$.
Conditioning on the state of the network at time $t$,

$p_{j}(t+\Delta t)=\sum_{x\in E}\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]\mathbb{P}^{\eta}[X_{t}=x]=\sum_{x\in
E}\bigl{(}\operatorname{Re}x_{j}+\operatorname{Im}x_{j}+1-\lvert
x_{j}\rvert\bigr{)}\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]\mathbb{P}^{\eta}[X_{t}=x].$

Because $x_{j}$ has a value in $\\{0,1,i\\}$, exactly one of
$\operatorname{Re}x_{j}$, $\operatorname{Im}x_{j}$ or $1-\lvert x_{j}\rvert$
is 1 while the other two are 0. Hence, we can now split the sum in three
parts: one for each possible state of $j$ at time $t$. Each conditional
probability $\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]$ can then be expressed
from transition rates according to the relations given in eq. (1):

$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=0]$
$\displaystyle=a_{j}(\eta,x)\Delta t+o(\Delta t),$
$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=1]$
$\displaystyle=1-\beta_{j}(\eta)\Delta t+o(\Delta t),$
$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=i]$
$\displaystyle=o(\Delta t).$

Using these relations, we obtain that as $\Delta t\downarrow 0$,

$\displaystyle p_{j}(t+\Delta t)$ $\displaystyle=\sum_{x\in
E}\operatorname{Re}x_{j}\bigl{(}1-\beta_{j}(\eta)\Delta
t\bigr{)}\mathbb{P}^{\eta}[X_{t}=x]+\sum_{x\in E}\bigl{(}1-\lvert
x_{j}\rvert\bigr{)}a_{j}(\eta,x)\Delta t\mathbb{P}^{\eta}[X_{t}=x]+o(\Delta
t).$ These sums are now expectations of functions of $X_{t}$, and we find that
$\displaystyle p_{j}(t+\Delta t)$
$\displaystyle=\bigl{(}1-\beta_{j}(\eta)\Delta t\bigr{)}p_{j}(t)+\Delta
t\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}a_{j}(\eta,X_{t})\big{]}+o(\Delta t).$

as $\Delta t\downarrow 0$. Rearranging and dividing through by $\Delta t$
yields

$\dot{p}_{j}(t)=-\beta_{j}(\eta)p_{j}(t)+\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}a_{j}(\eta,X_{t})\big{]}.$

The derivatives of $r_{j}$ and $q_{j}$ are obtained in the same way. In these
two cases, the transition probabilities $\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]$ and
$\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=0\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]$ are expressed from other
transition rates, and we obtain

$\dot{r}_{j}(t)=-\gamma_{j}(\eta)r_{j}(t)+\beta_{j}(\eta)p_{j}(t)\qquad\text{and}\qquad\dot{q}_{j}(t)=-\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}a_{j}(\eta,X_{t})\big{]}+\gamma_{j}(\eta)r_{j}(t).$

### A.2 The nonautonomous macroscopic first-order evolution equation

Now, we use the expressions of the derivatives of $p_{j}$, $q_{j}$ and $r_{j}$
obtained above to find the derivatives of the expected fractions of
populations

$\E{A}[J](t)\vcentcolon=\mathbb{E}^{\eta}[A_{t}^{J}],\qquad\E{R}[J](t)\vcentcolon=\mathbb{E}^{\eta}[R_{t}^{J}]\qquad\text{and}\qquad\E{S}[J](t)\vcentcolon=\mathbb{E}^{\eta}[S_{t}^{J}],$

where $A_{t}^{J}$, $R_{t}^{J}$ and $S_{t}^{J}$ are respectively the fractions
of active, refractory and sensitive neurons in population $J$ at time $t$.

To obtain an expression for the derivative of $\E{A}[J]$, we use the idea
described in the main text: we average the derivatives $\dot{p}_{j}$ over $J$
using the linearity of expectations and derivatives. From the expression of
$\dot{p}_{j}$, this first results in

$\dE{A}[J](t)=\mathbb{E}^{\eta}\bigg{[}-\frac{1}{\lvert J\rvert}\sum_{j\in
J}\beta_{j}(\eta)\operatorname{Re}X_{t}^{j}+\frac{1}{\lvert J\rvert}\sum_{j\in
J}a_{j}(\eta,X_{t})\bigl{(}1-\lvert X_{t}^{j}\rvert\bigr{)}\Bigr{)}\bigg{]}.$

Introducing for $\xi\in\\{0,1,i\\}$ the subpopulations
$J_{t}^{\xi}\vcentcolon=\\{j\in J:X_{t}^{j}=\xi\\}$, the last equality can be
written as

$\dE{A}[J](t)=\mathbb{E}^{\eta}\bigg{[}-\frac{A_{t}^{J}}{\lvert
J_{t}^{1}\rvert}\sum_{j\in J_{t}^{1}}\beta_{j}(\eta)+\frac{S_{t}^{J}}{\lvert
J_{t}^{0}\rvert}\sum_{j\in J_{t}^{0}}a_{j}(\eta,X_{t})\bigg{]},$ (21)

where we used the identities $\lvert J_{t}^{1}\rvert=\lvert J\rvert A_{t}^{J}$
and $\lvert J_{t}^{0}\rvert=\lvert J\rvert S_{t}^{J}$. Now, assuming that the
number of neurons in $J$ is large, we expect the number of neurons in $J$ that
are active at time $t$ to be large as well. Since the random variables
$\beta_{j}$ for $j\in J_{t}^{1}$ are independent and identically distributed,
the law of large numbers motivates the approximation

$\frac{1}{\lvert J_{t}^{1}\rvert}\sum_{j\in
J_{t}^{1}}\beta_{j}(\eta)\approx\mathbb{E}_{\mu}[\beta_{j}]=\vcentcolon\beta_{J},$
(22)

where $\mathbb{E}_{\mu}$ denotes the expectation on $(H,\mathscr{H},\mu)$. The
other term can be handled in a similar way. First, we approximate the input

$\sum_{k=1}^{N}W_{jk}(\eta)\operatorname{Re}X_{t}^{k}+Q_{J}=\sum_{K\in\mathscr{P}}\sum_{k\in
K_{t}^{1}}W_{jk}(\eta)+Q_{J}\approx\sum_{K\in\mathscr{P}}\lvert K\rvert
A_{t}^{K}\mathbb{E}_{\mu}[W_{jk}]+Q_{J}.$

To simplify notation, we define the input in population $J$ at time $t$

$B_{t}^{J}\vcentcolon=\sum_{K\in\mathscr{P}}c_{JK}A_{t}^{K}+Q_{J}\quad\text{with}\quad
c_{JK}\vcentcolon=\lvert K\rvert\mathbb{E}_{\mu}[W_{jk}]$

for $j\in J$ and $k\in K$. This leads to approximate, for $j\in J$,

$a_{j}(\eta,X_{t})\approx\alpha_{j}(\eta){}_{\\{B_{t}^{J}>\theta_{j}\\}}(\eta),$

and the law of large numbers now motivates the approximation

$\frac{1}{\lvert J_{t}^{0}\rvert}\sum_{j\in
J_{t}^{0}}a_{j}(\eta,X_{t})\approx\mathbb{E}_{\mu}[\alpha_{j}{}_{\\{B_{t}^{J}>\theta_{j}\\}}]=\alpha_{J}F_{\theta_{J}}(B_{t}^{J}),$
(23)

where $\alpha_{J}\vcentcolon=\mathbb{E}_{\mu}[\alpha_{j}]$ and
$F_{\theta_{J}}$ denotes the cumulative distribution function of $\theta_{j}$
for $j\in J$.

Using the approximations from eq. (22) and eq. (23) in the expression of
$\dE{A}[J]$ given in eq. (21) yields

$\dE{A}[J](t)\approx-\beta_{J}\E{A}[J](t)+\alpha_{J}\mathbb{E}^{\eta}[S_{t}^{J}F_{\theta_{J}}(B_{t}^{J})].$

The same method can be applied to find approximate expressions for the
derivatives of $\E{R}[J]$ and $\E{S}[J]$. This leads us to model the
macroscopic dynamics of the network by the differential equations given in eq.
(5) in the main text.

### A.3 The nonautonomous microscopic second-order evolution equation

As a first step towards finding the derivatives of covariances given in eq.
(7), we now find expressions for the derivatives of the probabilities

$\displaystyle\SwapAboveDisplaySkip p_{jk}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=X_{t}^{k}=1]=\mathbb{E}^{\eta}[\operatorname{Re}X_{t}^{j}\operatorname{Re}X_{t}^{k}],$
$\displaystyle r_{jk}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=X_{t}^{k}=i]=\mathbb{E}^{\eta}[\operatorname{Im}X_{t}^{j}\operatorname{Im}X_{t}^{k}],$
$\displaystyle\rho_{jk}(t)$
$\displaystyle\vcentcolon=\mathbb{P}^{\eta}[X_{t}^{j}=1,X_{t}^{k}=i]=\mathbb{E}^{\eta}[\operatorname{Re}X_{t}^{j}\operatorname{Im}X_{t}^{k}].$

The idea is exactly the same as in the case of the probabilities $p_{j}$,
$r_{j}$ and $q_{j}$. We will give the details of the calculations only for the
case of $\rho_{jk}$, as it illustrates all of the relevant ideas.

Let $\Delta t>0$. Then

$\rho_{jk}(t+\Delta t)=\sum_{x\in
E}\bigl{(}\operatorname{Re}x_{j}+\operatorname{Im}x_{j}+1-\lvert
x_{j}\rvert\bigr{)}\bigl{(}\operatorname{Re}x_{k}+\operatorname{Im}x_{k}+1-\lvert
x_{k}\rvert\bigr{)}\mathbb{P}^{\eta}[X_{t+\Delta t}^{j}=1,X_{t+\Delta
t}^{k}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]\mathbb{P}^{\eta}[X_{t}=x].$

The sum now gets split into nine parts, one for each of the possible states
for the pair $(j,k)$. Again, the conditional probability
$\mathbb{P}^{\eta}[X_{t+\Delta t}^{j}=1,X_{t+\Delta
t}^{k}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]$ can be expressed in
terms of transition rates. These expressions are all derived from the general
statement that for any $x,y\in E$, as $\Delta t\downarrow 0$,

$\mathbb{P}^{\eta}[X_{t+\Delta
t}=y\nonscript\,|\nonscript\,\mathopen{}X_{t}=x]=\delta_{xy}+m^{\eta}(x,y)\Delta
t+o(\Delta t),$

where $m^{\eta}(x,y)$ is the $(x,y)$ entry of the generator $M^{\eta}$, as in
the main text. It follows from this statement that as $\Delta t\downarrow 0$,

$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta t}^{j}=1,X_{t+\Delta
t}^{k}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=1,X_{t}^{k}=i]$
$\displaystyle=1-\bigl{(}\beta_{j}(\eta)+\gamma_{k}(\eta)\bigr{)}\Delta
t+o(\Delta t),$ $\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta
t}^{j}=1,X_{t+\Delta
t}^{k}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=0,X_{t}^{k}=i]$
$\displaystyle=a_{j}(\eta,x)\Delta t+o(\Delta t),$
$\displaystyle\mathbb{P}^{\eta}[X_{t+\Delta t}^{j}=1,X_{t+\Delta
t}^{k}=i\nonscript\,|\nonscript\,\mathopen{}X_{t}=x,X_{t}^{j}=1,X_{t}^{k}=1]$
$\displaystyle=\beta_{k}(\eta)\Delta t+o(\Delta t),$

while the transition probabilities from states $x$ where $j$ and $k$ have
other states are all $o(\Delta t)$. Using these relations, we see that as
$\Delta t\downarrow 0$,

$\displaystyle\rho_{jk}(t+\Delta t)$ $\displaystyle=\sum_{x\in
E}\operatorname{Re}x_{j}\operatorname{Im}x_{k}\bigl{(}1-\bigl{(}\beta_{j}(\eta)+\gamma_{k}(\eta)\bigr{)}\Delta
t\bigr{)}\mathbb{P}^{\eta}[X_{t}=x]$ $\displaystyle\hskip
62.59605pt+\sum_{x\in E}\bigl{(}1-\lvert
x_{j}\rvert\bigr{)}\operatorname{Im}x_{k}\,a_{j}(\eta,x)\Delta
t\mathbb{P}^{\eta}[X_{t}=x]+\sum_{x\in
E}\operatorname{Re}x_{j}\operatorname{Re}x_{k}\,\beta_{k}(\eta)\Delta
t\mathbb{P}^{\eta}[X_{t}=x]+o(\Delta t)$
$\displaystyle=\bigl{(}1-\bigl{(}\beta_{j}(\eta)+\gamma_{k}(\eta)\Delta
t\bigr{)}\rho_{jk}(t)+\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\operatorname{Im}X_{t}^{k}a_{j}(\eta,X_{t})\big{]}\Delta
t+\beta_{k}(\eta)p_{jk}(t)\Delta t+o(\Delta t).$

Rearranging and dividing through by $\Delta t$ finally yields

$\displaystyle\dot{\rho}_{jk}(t)$
$\displaystyle=-\bigl{(}\beta_{j}(\eta)+\gamma_{k}(\eta)\bigr{)}\rho_{jk}(t)+\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\operatorname{Im}X_{t}^{k}a_{j}(\eta,X_{t})\big{]}+\beta_{k}(\eta)p_{jk}(t).$
(24a) The same method allows to find $\displaystyle\dot{p}_{jk}(t)$
$\displaystyle=-\bigl{(}\beta_{j}(\eta)+\beta_{k}(\eta)\bigr{)}p_{jk}(t)+\mathbb{E}^{\eta}\big{[}\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\operatorname{Re}X_{t}^{k}a_{j}(\eta,X_{t})\big{]}+\mathbb{E}^{\eta}\big{[}\operatorname{Re}X_{t}^{j}\bigl{(}1-\lvert
X_{t}^{k}\rvert\bigr{)}a_{k}(\eta,X_{t})\big{]},$ (24b)
$\displaystyle\dot{r}_{jk}(t)$
$\displaystyle=-\bigl{(}\gamma_{j}(\eta)+\gamma_{k}(\eta)\bigr{)}r_{jk}(t)+\beta_{j}(\eta)\rho_{jk}(t)+\beta_{k}(\eta)\rho_{kj}(t).$
(24c)

We could find similar equations to describe the evolution of similar
probabilities with $j$ or $k$ being sensitive, but these can always be
expressed as functions of $p_{jk}$, $r_{jk}$, $\rho_{jk}$ and $\rho_{kj}$
since each entry $X_{t}^{j}$ is always in $\\{0,1,i\\}$.

### A.4 The nonautonomous macroscopic second-order evolution equation

We now find the derivatives of covariances between active and refractory
fractions of populations of the network. We first find the derivatives of the
expectations

$\EE{AA}[JK](t)\vcentcolon=\mathbb{E}^{\eta}[A_{t}^{J}A_{t}^{K}],\qquad\EE{RR}[JK](t)\vcentcolon=\mathbb{E}^{\eta}[R_{t}^{J}R_{t}^{K}]\qquad\text{and}\qquad\EE{AR}[JK](t)\vcentcolon=\mathbb{E}^{\eta}[A_{t}^{J}R_{t}^{K}]$

from the differential equations eq. (24), using the same strategy as for first
moments. As in the microscopic case we only detail the case of $\EE{AR}[JK]$,
as the other two are similar.

By linearity of expectations and derivatives,

$\displaystyle\dEE{AR}[JK](t)$
$\displaystyle=\mathbb{E}^{\eta}\bigg{[}\frac{1}{\lvert J\rvert\lvert
K\rvert}\sum_{j\in J,k\in
K}\Bigl{(}-\bigl{(}\beta_{j}(\eta)+\gamma_{k}(\eta)\bigr{)}\operatorname{Re}X_{t}^{j}\operatorname{Im}X_{t}^{k}+\bigl{(}1-\lvert
X_{t}^{j}\rvert\bigr{)}\operatorname{Im}X_{t}^{k}a_{j}(\eta,X_{t})+\beta_{k}(\eta)\operatorname{Re}X_{t}^{j}\operatorname{Re}X_{t}^{k}\Bigr{)}\bigg{]}$
$\displaystyle=\mathbb{E}^{\eta}\bigg{[}-\frac{A_{t}^{J}R_{t}^{K}}{\lvert
J_{t}^{1}\rvert}\sum_{j\in
J_{t}^{1}}\beta_{j}(\eta)-\frac{A_{t}^{J}R_{t}^{K}}{\lvert
K_{t}^{i}\rvert}\sum_{k\in
K_{t}^{i}}\gamma_{k}(\eta)+\frac{S_{t}^{J}R_{t}^{K}}{\lvert
J_{t}^{0}\rvert}\sum_{j\in
J_{t}^{0}}a_{j}(\eta,X_{t})+\frac{A_{t}^{J}A_{t}^{K}}{\lvert
K_{t}^{1}\rvert}\sum_{k\in K_{t}^{1}}\beta_{k}(\eta)\bigg{]}.$

In the same way as in the case of first-order moments, the law of large
numbers motivates the approximations

$\frac{1}{\lvert J_{t}^{1}\rvert}\sum_{j\in
J_{t}^{1}}\beta_{j}(\eta)\approx\beta_{J},\qquad\frac{1}{\lvert
K_{t}^{i}\rvert}\sum_{k\in
K_{t}^{i}}\gamma_{k}(\eta)\approx\gamma_{K}\qquad\text{and}\qquad\frac{1}{\lvert
J_{t}^{0}\rvert}\sum_{j\in
J_{t}^{0}}a_{j}(\eta,X_{t})\approx\alpha_{J}F_{\theta_{J}}(B_{t}^{J}),$

and we find the approximate equation

$\displaystyle\dEE{AR}[JK](t)$
$\displaystyle=-(\beta_{J}+\gamma_{K})\EE{AR}[JK](t)+\alpha_{J}\mathbb{E}^{\eta}[S_{t}^{J}R_{t}^{K}F_{\theta_{J}}(B_{t}^{J})]+\beta_{K}\EE{AA}[JK](t).$
(25a) The same method leads to $\displaystyle\dEE{AA}[JK](t)$
$\displaystyle=-(\beta_{J}+\beta_{K})\EE{AA}[JK](t)+\alpha_{J}\mathbb{E}^{\eta}[S_{t}^{J}A_{t}^{K}F_{\theta_{J}}(B_{t}^{J})]+\alpha_{K}\mathbb{E}^{\eta}[A_{t}^{J}S_{t}^{K}F_{\theta_{K}}(B_{t}^{K})],$
(25b) $\displaystyle\dEE{RR}[JK](t)$
$\displaystyle=-(\gamma_{J}+\gamma_{K})\EE{RR}[JK](t)+\beta_{J}\EE{AR}[JK](t)+\beta_{K}\EE{AR}[KJ](t).$
(25c)

By definition of the covariance, $\C{AR}[JK]=\EE{AR}[JK]-\E{A}[J]\E{R}[K]$,
and similar relations hold for other combinations of fractions of populations.
Therefore, the expressions of the derivatives of the covariances $\C{AA}[JK]$,
$\C{RR}[JK]$ and $\C{AR}[JK]$ given in eq. (7) follow from eq. (25) along with
the derivatives of $\E{A}[J]$ and $\E{R}[J]$ given in eq. (5).

### A.5 Limits of the approximations

We start by remarking that, even though we motivate the approximations of the
averages of transition rates (e.g. in eq. (22) and eq. (23)) by the law of
large numbers, it would not be obvious to make the argument fully rigorous by
taking a limit where the sizes of all populations grow infinitely large. The
reason for this is that the sets that are needed to be large, which are the
subpopulations of active, sensitive and refractory neurons in each population,
are themselves random on $(\Omega,\mathscr{F},\mathbb{P}^{\eta})$, but the
probability measure $\mathbb{P}^{\eta}$ is a function on
$(H,\mathscr{H},\mu)$. Therefore, the setting cannot be easily translated to a
simple sequence of independent and identically distributed random variables,
and we settle for understanding eq. (22) and eq. (23) as approximations.

We expect the approximations

$\frac{1}{\lvert J_{t}^{1}\rvert}\sum_{j\in
J_{t}^{1}}\beta_{j}(\eta)\approx\beta_{J}\qquad\text{and}\qquad\frac{1}{\lvert
J_{t}^{i}\rvert}\sum_{j\in J_{t}^{i}}\gamma_{j}(\eta)\approx\gamma_{J}$

to be good approximations for large populations regardless of the specific
distributions of $\beta_{j}$ and $\gamma_{j}$ in $J$, as long as these remain
independent and identically distributed in $J$. Indeed, we do not expect the
specific values of $\beta_{j}$ and $\gamma_{j}$ to have an impact on the
probabilities for $j$ to be active and refractory at time $t$, so we expect
the empirical distributions of $\beta_{j}$ and $\gamma_{j}$ over $J_{t}^{1}$
and $J_{t}^{i}$ to approximate reasonably well the actual distributions when
$J$ is large.

However, the case of the corresponding approximations for the activation
rates,

$\frac{1}{\lvert J_{t}^{0}\rvert}\sum_{j\in
J_{t}^{0}}\alpha_{j}(\eta){}_{T_{j}(X_{t})}(\eta)\approx\frac{1}{\lvert
J_{t}^{0}\rvert}\sum_{j\in
J_{t}^{0}}\alpha_{j}(\eta){}_{\\{B_{t}^{J}>\theta_{j}\\}}(\eta)\approx\alpha_{J}F_{\theta_{J}}(B_{t}^{J}),$
(26)

requires more care. Assuming that every neuron of $J$ is given the same input
$B_{t}^{J}$, the last approximation should hold only when the distribution of
$\theta_{j}$ in $J$ is not too spread out. Indeed, only the neurons in $J$
whose thresholds are lower than their input have a nonzero probability to
activate, so it is possible for the specific value of the threshold of a
neuron to have an effect on the probability that it is sensitive at a given
time. This could lead the empirical distribution of the thresholds in the
subpopulation $J_{t}^{0}$ to be biased towards higher thresholds.
Nevertheless, since the activation probabilities only depend on the sign of
the difference between the input and the threshold, we do not expect this
effect to be important when the distribution of the thresholds is peaked
enough. Similar arguments also imply that the first approximation in eq. (26)
should only hold when the distributions of the weights $W_{jk}$ are peaked
enough.

##  Appendix B Proof of Theorem 1 and its corollaries

In this section, we prove Theorem 1 and its corollaries. For convenience, we
recall the general setup.

We study here a function $F_{\theta_{J}}\colon\mathbb{R}\to[0,1]$, which is
the cumulative distribution function of the thresholds in a population $J$ of
the network. We assume that this distribution is unimodal and symmetric with
mean $\theta_{J}$. Then, we define the function
$g_{J}\colon\mathbb{R}\times[0,\infty)$ by setting

$g_{J}(b,v)\vcentcolon=\frac{v}{2(\theta_{J}-b)}\frac{F_{\theta_{J}}^{\prime\prime}(b)}{F_{\theta_{J}}^{\prime}(b)}$
(27)

for $b\neq\theta_{J}$, assuming $F_{\theta_{J}}^{\prime}>0$. Since the
distribution is unimodal and symmetric around $\theta_{J}$, it must be that
$F_{\theta_{J}}^{\prime\prime}(\theta_{J})=0$, and $g_{J}$ can be continuously
extended to $b=\theta_{J}$ by definition of the derivative of
$F_{\theta_{J}}^{\prime\prime}$, as long as it exists. Finally, we define
$G_{J}\colon\mathbb{R}\times[0,\infty)$ as

$G_{J}(b,v)\vcentcolon=F_{\theta_{J}}\Bigl{(}\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\Bigr{)}.$
(28)

We now recall Theorem 1 and prove it.

######  Theorem 1.

Suppose that the thresholds in population $J$ follow a unimodal and symmetric
distribution with mean $\theta_{J}$ and cumulative distribution function
$F_{\theta_{J}}$. Let $g_{J}$ and $G_{J}$ be defined by eq. (27) and eq. (28)
respectively. Suppose that

  1. 1.

$F_{\theta_{J}}$ is $\mathscr{C}^{4}$ on $\mathbb{R}$;

  2. 2.

$F_{\theta_{J}}^{\prime}>0$;

  3. 3.

$g_{J}(\cdot,1)$ is bounded on $\mathbb{R}$;

  4. 4.

$\forall b\in\mathbb{R}$, $g_{J}(b,1)+(\theta_{J}-b)\partial_{1}g_{J}(b,1)\geq
0$.

Then $G_{J}$ satisfies the following conditions.

  1. i.

$G_{J}$ is $\mathscr{C}^{1}$ on $\mathbb{R}\times[0,\infty)$;

  2. ii.

$G_{J}(\cdot,0)=F_{\theta_{J}}$ and $\forall
b\in\mathbb{R},G_{J}(b,v)\to\nicefrac{{1}}{{2}}$ as $v\to\infty$;

  3. iii.

$\forall v\geq 0,G_{J}(b,v)\to 0$ as $b\to-\infty$ and $G_{J}(b,v)\to 1$ as
$b\to\infty$;

  4. iv.

$\forall v\geq 0$, $G_{J}(\cdot,v)$ is increasing;

  5. v.

$\partial_{2}G_{J}(\cdot,0)=\frac{1}{2}F_{\theta_{J}}^{\prime\prime}$.

###### Proof.

We start by proving the second and third properties. First, note that for any
$b\in\mathbb{R}$, $g_{J}(b,0)=0$ so $G_{J}(b,0)=F_{\theta_{J}}(b)$. Then, as
$v\to\infty$, $g_{J}(b,v)=vg_{J}(b,1)\to\infty$, which implies that

$\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\to\theta_{J}\qquad\text{so
that}\qquad G_{J}(b,v)\to F_{\theta_{J}}(\theta_{J})=\frac{1}{2},$

proving ii. Similarly, the fact that $g_{J}(\cdot,1)$ is bounded implies that
$g_{J}$ is bounded with respect to its first argument, so

$\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\to\pm\infty$

as $b\to\pm\infty$, and property iii follows by properties of a cumulative
distribution function.

To prove the other three properties, we compute the derivatives of $G_{J}$. A
simple calculation shows that

$\partial_{1}g_{J}(b,v)=\frac{v}{2(\theta_{J}-b)}\biggl{(}\frac{F_{\theta_{J}}^{\prime\prime\prime}(b)}{F_{\theta_{J}}^{\prime}(b)}+\frac{F_{\theta_{J}}^{\prime\prime}(b)}{F_{\theta_{J}}^{\prime}(b)(\theta_{J}-b)}-\Bigl{(}\frac{F_{\theta_{J}}^{\prime\prime}(b)}{F_{\theta_{J}}^{\prime}(b)}\Bigr{)}^{2}\biggr{)}$

for $b\neq\theta_{J}$, and applying l’Hospital’s rule shows that as
$b\to\theta_{J}$,
$\partial_{1}g_{J}(b,v)\to-\frac{F_{\theta_{J}}^{\prime\prime\prime\prime}(b)}{2F_{\theta_{J}}^{\prime}(\theta_{J})}$.
This is zero since $F_{\theta_{J}}^{\prime\prime\prime\prime}(\theta_{J})=0$,
as the distribution is unimodal and symmetric. It is easy to verify by
applying l’Hospital’s rule twice that $\partial_{1}g_{J}(\theta_{J},v)=0$, so
that $\partial_{1}g_{J}$ exists and is continuous on
$\mathbb{R}\times[0,\infty)$. As it is clear that $\partial_{2}g_{J}$ exists
and is continuous on $\mathbb{R}\times[0,\infty)$, it follows that $g_{J}$ is
$\mathscr{C}^{1}$. This directly implies that $G_{J}$ is $\mathscr{C}^{1}$ as
well.

Now, a direct computation leads to

$\partial_{1}G_{J}(b,v)=F_{\theta_{J}}^{\prime}\Bigl{(}\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\Bigr{)}\frac{1+g_{J}(b,v)+(\theta_{J}-b)\partial_{1}g_{J}(b,v)}{\bigl{(}1+g_{J}(b,v)\bigr{)}^{2}}.$

This is always positive by assumption 4, and property iv follows. We record
here for future computations that the above implies that

$\partial_{1}G_{J}(b,0)=F_{\theta_{J}}^{\prime}(b).$ (29)

Finally,

$\partial_{2}G_{J}(b,v)=F_{\theta_{J}}^{\prime}\Bigl{(}\frac{b+\theta_{J}g_{J}(b,v)}{1+g_{J}(b,v)}\Bigr{)}\frac{(\theta_{J}-b)\partial_{2}g_{J}(b,v)}{\bigl{(}1+g_{J}(b,v)\bigr{)}^{2}},$

and property v follows by evaluating at $v=0$. ∎

We now prove the corollaries of the Theorem.

######  Corollary 1.

Suppose that all assumptions of Theorem 1 hold. Then
$f\colon(0,\infty)\times\mathbb{R}^{3}\to(0,\infty)$ defined by

$f(s,b,c,v)\vcentcolon=sG_{J}\Bigl{(}b+\frac{c}{s},v\Bigr{)}$

satisfies the following conditions for any $s>0$, $b,c\in\mathbb{R}$ and
$v\geq 0$.

  1. i.

$f$ is $\mathscr{C}^{1}$ on $(0,\infty)\times\mathbb{R}^{3}$;

  2. ii.

$f(s,b,c,v)\gtreqless sG_{J}(b,v)$ when $c\gtreqless 0$, and in particular
$f(s,b,0,v)=sG_{J}(b,v)$ and $f(s,b,0,0)=sF_{\theta_{J}}(b)$;

  3. iii.

$f(s,b,c,v)\to 0$ as $b\to-\infty$ and $f(s,b,c,v)\to s$ as $b\to\infty$;

  4. iv.

$\partial_{3}f(s,b,0,0)=F_{\theta_{J}}^{\prime}(b)$;

  5. v.

$\partial_{4}f(s,b,0,0)=\frac{1}{2}sF_{\theta_{J}}^{\prime\prime}(b)$.

###### Proof.

Properties i, ii and iii follow respectively from properties i, iv and iii of
$G_{J}$ given in Theorem 1. Then, properties iv and v follow by evaluating the
derivatives using eq. (29) and property v of $G_{J}$. ∎

######  Corollary 2.

Suppose that all assumptions of Theorem 1 hold. Then
$f\colon(0,\infty)^{2}\times\mathbb{R}^{5}\to\mathbb{R}$ defined by

$f(x,s,b,c_{1},c_{2},c_{3},v)\vcentcolon=(xs+c_{1})G_{J}\Bigl{(}b+\frac{c_{2}}{x}+\frac{c_{3}}{s},v\Bigr{)}$

satisfies the following conditions for any $x,s>0$, $b\in\mathbb{R}$ and
$v\geq 0$.

  1. i.

$f$ is $\mathscr{C}^{1}$ on $(0,\infty)^{2}\times\mathbb{R}^{5}$;

  2. ii.

$f(x,s,b,0,0,0,v)=xsG_{J}(b,v)$, and in particular
$f(x,s,b,0,0,0,0)=xsF_{\theta_{J}}(b)$;

  3. iii.

$\partial_{4}f(x,s,b,0,0,0,0)=F_{\theta_{J}}(b)$;

  4. iv.

$\partial_{5}f(x,s,b,0,0,0,0)=sF_{\theta_{J}}^{\prime}(b)$;

  5. v.

$\partial_{6}f(x,s,b,0,0,0,0)=xF_{\theta_{J}}^{\prime}(b)$;

  6. vi.

$\partial_{7}f(x,s,b,0,0,0,0)=\frac{1}{2}xsF_{\theta_{J}}^{\prime\prime}(b)$.

###### Proof.

Property i follows from property i of $G_{J}$ in Theorem 1. Then, property ii
follows by direct evaluation, and properties iii to vi follow by evaluating
the derivatives using eq. (29) and properties ii and v of $G_{J}$. ∎

### B.1 Verification of the assumptions for a normal distribution

Suppose that the thresholds in population $J$ follow a normal distribution
with mean $\theta_{J}$ and variance $\sigma_{\theta_{J}}^{2}$, which has
density

$F_{\theta_{J}}^{\prime}(b)=\frac{1}{\sigma_{\theta_{J}}\sqrt{2\pi}}\exp\Bigl{(}-\frac{(b-\theta_{J})^{2}}{2\sigma_{\theta_{J}}^{2}}\Bigr{)}.$

This distribution is indeed unimodal and symmetric. Moreover, its cumulative
distribution function $F_{\theta_{J}}$ is smooth, and the density
$F_{\theta_{J}}^{\prime}$ is positive, so assumptions 1 and 2 are satisfied.

To verify the other two assumptions, we compute $g_{J}$ for this distribution.

$F_{\theta_{J}}^{\prime\prime}(b)=\frac{\theta_{J}-b}{\sigma_{\theta_{J}}^{3}\sqrt{2\pi}}\exp\Bigl{(}-\frac{(b-\theta_{J})^{2}}{2\sigma_{\theta_{J}}^{2}}\Bigr{)}=\frac{\theta_{J}-b}{\sigma_{\theta_{J}}^{2}}F_{\theta_{J}}^{\prime}(b),\qquad\text{so}\qquad
g_{J}(b,v)=\frac{v}{2\sigma_{\theta_{J}}^{2}}.$

Then $g_{J}(\cdot,1)\equiv\nicefrac{{1}}{{2\sigma_{\theta_{J}}^{2}}}$ is
bounded so assumptions 3 is satisfied, and it is positive so 4 is directly
satisfied since $g_{J}$ does not depend on its first argument.

Therefore, the conclusions of Theorem 1 hold if the thresholds in population
$J$ follow a normal distribution.

### B.2 Verification of the assumptions for a logistic distribution

Suppose that the thresholds in population $J$ follow a logistic distribution
with mean $\theta_{J}$ and scaling factor $s_{\theta_{J}}$, which has
cumulative distribution function

$F_{\theta_{J}}(b)=\sigma\Bigl{(}\frac{b-\theta_{J}}{s_{\theta_{J}}}\Bigr{)}\qquad\text{where}\qquad\sigma(x)\vcentcolon=\frac{1}{1+e^{-x}}.$

This distribution is unimodal and symmetric, and $F_{\theta_{J}}$ is smooth so
assumption 1 is satisfied.

The logistic function $\sigma$ has the property that
$\sigma^{\prime}=\sigma(1-\sigma)$. This property allows to compute

$F_{\theta_{J}}^{\prime}=\frac{1}{s_{\theta_{J}}}F_{\theta_{J}}(1-F_{\theta_{J}})\qquad\text{and}\qquad
F_{\theta_{J}}^{\prime\prime}=\frac{1}{s_{\theta_{J}}}F_{\theta_{J}}^{\prime}(1-2F_{\theta_{J}}).$

Since the range of $F_{\theta_{J}}$ is $(0,1)$, it follows that
$F_{\theta_{J}}^{\prime}>0$ and assumption 2 is satisfied. Then, for
$b\neq\theta_{J}$,

$g_{J}(b,v)=\frac{v}{2s_{\theta_{J}}}\frac{1-2F_{\theta_{J}}(b)}{\theta_{J}-b},$

and $g_{J}$ is extended to a continuous function on
$\mathbb{R}\times[0,\infty)$ using l’Hospital’s rule. Since $g_{J}(\cdot,1)$
is continuous on $\mathbb{R}$, in order to prove that it is bounded on
$\mathbb{R}$ it suffices to prove that it is bounded outside of the compact
interval $[\theta_{J}-1,\theta_{J}+1]$, and this is easy to see: if $\lvert
b-\theta_{J}\rvert>1$, then $\lvert
g_{J}(b,1)\rvert<\nicefrac{{1}}{{2s_{\theta_{J}}}}$ since $F_{\theta_{J}}$
takes its values in $(0,1)$. This shows that assumption 3 is satisfied.

To verify the last assumption, notice that for $b\neq\theta_{J}$ and any
$v\geq 0$,

$(\theta_{J}-b)\partial_{1}g_{J}(b,v)=\frac{v}{2s_{\theta_{J}}}\Bigl{(}\frac{1-2F_{\theta_{J}}(b)}{\theta_{J}-b}-2F_{\theta_{J}}^{\prime}(b)\Bigr{)}=g_{J}(b,v)\Bigl{(}1-\frac{2(\theta_{J}-b)F_{\theta_{J}}^{\prime}(b)}{1-2F_{\theta_{J}}(b)}\Bigr{)}.$

Then, with $x\vcentcolon=\nicefrac{{(b-\theta_{J})}}{{s_{\theta_{J}}}}$,

$1-2F_{\theta_{J}}(b)=-\frac{1-e^{-x}}{1+e^{-x}}\qquad\text{and}\qquad
F_{\theta_{J}}^{\prime}(b)=\frac{1}{s_{\theta_{J}}}\frac{e^{-x}}{(1+e^{-x})^{2}},$

so

$(\theta_{J}-b)\partial_{1}g_{J}(b,v)=g_{J}(b,v)\Bigl{(}1-\frac{2xe^{-x}}{1-e^{-2x}}\Bigr{)}=g_{J}(b,v)\Bigl{(}1-\frac{x}{\sinh
x}\Bigr{)}.$

Since $\lvert x\rvert\leq\lvert\sinh x\rvert$ for all $x\in\mathbb{R}$, it
follows that $(\theta_{J}-b)\partial_{1}g_{J}(b,v)\geq 0$ because $g_{J}$ is
nonnegative. Since this holds whenever $b\neq\theta_{J}$, it holds for
$b=\theta_{J}$ by continuity, and assumption 4 holds.

Therefore, the conclusions of Theorem 1 hold if the thresholds in population
$J$ follow a logistic distribution.