# Macroscopic description for networks of spiking neurons

Ernest Montbrió  Center for Brain and Cognition. Department of Information and
Communication Technologies, Universitat Pompeu Fabra, 08018 Barcelona, Spain
Diego Pazó  Instituto de Física de Cantabria (IFCA), CSIC-Universidad de
Cantabria, 39005 Santander, Spain  Alex Roxin  Centre de Recerca Matemàtica,
Campus de Bellaterra, 08193 Bellaterra (Barcelona), Spain.

###### Abstract

A major goal of neuroscience, statistical physics and nonlinear dynamics is to
understand how brain function arises from the collective dynamics of networks
of spiking neurons. This challenge has been chiefly addressed through large-
scale numerical simulations. Alternatively, researchers have formulated mean-
field theories to gain insight into macroscopic states of large neuronal
networks in terms of the collective firing activity of the neurons, or the
firing rate. However, these theories have not succeeded in establishing an
exact correspondence between the firing rate of the network and the underlying
microscopic state of the spiking neurons. This has largely constrained the
range of applicability of such macroscopic descriptions, particularly when
trying to describe neuronal synchronization. Here we provide the derivation of
a set of exact macroscopic equations for a network of spiking neurons. Our
results reveal that the spike generation mechanism of individual neurons
introduces an effective coupling between two biophysically relevant
macroscopic quantities, the firing rate and the mean membrane potential, which
together govern the evolution of the neuronal network. The resulting equations
exactly describe all possible macroscopic dynamical states of the network,
including states of synchronous spiking activity. Finally we show that the
firing rate description is related, via a conformal map, with a low-
dimensional description in terms of the Kuramoto order parameter, called Ott-
Antonsen theory. We anticipate our results will be an important tool in
investigating how large networks of spiking neurons self-organize in time to
process and encode information in the brain.

###### pacs:

87.19.lj 05.45.Xt 87.10.-Ed 05.65.+b

Processing and coding of information in the brain necessarily imply the
coordinated activity of large ensembles of neurons. Within sensory regions of
the cortex, many cells show similar responses to a given stimulus, indicating
a high degree of neuronal redundancy at the local level. This suggests that
information is encoded in the population response and hence can be captured
via macroscopic measures of the network activity Averbeck  _et al._ (2006).
Moreover, the collective behavior of large neuronal networks is particularly
relevant given that current brain measurement techniques, such as
electroencephalography (EEG) or functional magnetic resonance imaging (fMRI),
provide data which is necessarily averaged over the activity of a large number
of neurons.

The macroscopic dynamics of neuronal ensembles has been extensively studied
through computational models of large networks of recurrently coupled spiking
neurons, including Hodgkin-Huxley-type conductance-based neurons Hodgkin and
Huxley (1952) as well as simplified neuron models, see e.g. Brette  _et al._
(2007); Izhikevich and Edelman (2008); Izhikevich (2007). In parallel,
researchers have sought to develop statistical descriptions of neuronal
networks, mainly in terms of a macroscopic observable that measures the mean
rate at which neurons emit spikes, the _firing rate_ Wilson and Cowan (1972);
Amari (1974); Freeman (1975); Gerstner (1995); Fusi and Mattia (1998); Amit
and Brunel (1997); Brunel (2000); Dayan and Abbott (2001); Gerstner and
Kistler (2002); Shriki  _et al._ (2003); Coombes (2005); Roxin  _et al._
(2005); Deco  _et al._ (2008); Ermentrout and Terman (2010); Ostojic and
Brunel (2011); Roxin and Montbrió (2011); Schaffer  _et al._ (2013). These
descriptions, called firing-rate equations (FREs), have been proven to be
extremely useful in understanding general computational principles underlying
functions such as memory Hopfield (1984); Mongillo  _et al._ (2008), visual
processing Ben-Yishai  _et al._ (1995); Hansel and Sompolinsky (1998); Moreno-
Bote  _et al._ (2007), motor control Zhang (1996) or decision making Wong and
Wang (2006).

Despite these efforts, to date there is no exact theory linking the dynamics
of a large network of spiking neurons with that of the firing rate.
Specifically, current macroscopic descriptions do not offer a precise
correspondence between the microscopic dynamics of the individual neurons,
e.g. individual membrane potentials, and the firing rate dynamics of the
neuronal network.

Indeed, FREs are generally derived through heuristic arguments which rely on
the underlying spiking activity of the neurons being asynchronous and hence
uncorrelated. As such, firing rate descriptions are not sufficient to describe
network states involving some degree of spike synchronization. Synchronization
is, however, an ubiquitous phenomenon in the brain, and its potential role in
neuronal computation is the subject of intense research Gerstner and Kistler
(2002); Abeles (1991); Varela  _et al._ (2001); Engel and Singer (2001); Ward
(2003); Fries (2005); Buzsaki (2006); Fell and Axmacher (2011). Hence, the
lack of firing rate descriptions for synchronous states limits the range of
applicability of mean-field theories to investigate neuronal dynamics.

Here we propose a method to derive the FREs for networks of heterogeneous,
all-to-all coupled quadratic integrate-and-fire (QIF) neurons, which is exact
in the thermodynamic limit, i.e. for large numbers of neurons. We consider an
ansatz for the distribution of the neurons’ membrane potentials that we
denominate the _Lorentzian ansatz_ (LA). The LA solves the corresponding
continuity equation exactly, making the system amenable to theoretical
analysis. Specifically, for particular distributions of the heterogeneity, the
LA yields a system of two ordinary differential equations for the firing rate
and mean membrane potential of the neuronal population. These equations fully
describe the macroscopic states of the network —including synchronized
states—, and represent the first example of an exact firing-rate description
for a network of recurrently connected spiking neurons. We finally show how
the LA transforms, via a conformal mapping, into the so-called Ott-Antonsen
ansatz (OA) that is extensively used to investigate the low-dimensional
dynamics of large populations of phase oscillators in terms of the Kuramoto
order parameter Ott and Antonsen (2008).

##  I Model description

Hodgkin-Huxley-type neuronal models can be broadly classified into two
classes, according to the nature of their transition to spiking in response to
an injected current Hodgkin (1948); Rinzel and Ermentrout (1998). Neuronal
models with so-called _Class I_ excitability, generate action potentials with
arbitrarily low frequency, depending on the strength of the applied current.
This occurs when a resting state disappears through a saddle-node bifurcation.
In contrast, in neurons with _Class II_ excitability the action potentials are
generated with a finite frequency. This occurs when the resting state loses
stability via a Hopf bifurcation. The QIF neuron is the canonical model for
_Class I_ neurons, and thus generically describes their dynamics near the
spiking threshold Ermentrout and Kopell (1986); Latham  _et al._ (2000);
Izhikevich (2007). Our aim here is to derive the FREs corresponding to a
heterogeneous all-to-all coupled population of $N$ QIF neurons. The
correspondence is exact in the thermodynamic limit, i.e. when $N\to\infty$
(this convergence has been recently studied in Buice and Chow (2013)).

The microscopic state of the population of QIF neurons is given by the
membrane potentials $\\{V_{j}\\}_{j=1,\ldots,N}$, which obey the following
ordinary differential equations Izhikevich (2007):

$\dot{V}_{j}=V_{j}^{2}+I_{j},\quad\text{if }V_{j}\geq V_{p},\text{then
}V_{j}\leftarrow V_{r}.$ (1)

Here the overdot denotes the time derivative, and $I_{j}$ represents an input
current. Each time a neuron’s membrane potential $V_{j}$ reaches the peak
value $V_{p}$, the neuron emits a spike and its voltage is reset to the value
$V_{r}$. In our analysis we consider the limit $V_{p}=-V_{r}\to\infty$. This
resetting rule captures the spike reset as well as the refractoriness of the
neuron. Without loss of generality we have rescaled the time and the voltage
in (1) to absorb any coefficients which would have appeared in the first two
terms. The form for the input currents is:

$I_{j}=\eta_{j}+Js(t)+I(t),$ (2)

where the external input has a heterogeneous, quenched component $\eta_{j}$ as
well as a common time-varying component $I(t)$, and the recurrent input is the
synaptic weight $J$ times the mean synaptic activation $s(t)$, which is
written:

$s(t)=\frac{1}{N}\sum_{j=1}^{N}\sum_{k\backslash
t_{j}^{k}<t}\int_{-\infty}^{t}dt^{\prime}a_{\tau}(t-t^{\prime})\delta(t^{\prime}-t_{j}^{k}).$
(3)

Here, $t_{j}^{k}$ is the time of the $k$th spike of $j$th neuron, $\delta(t)$
is the Dirac delta function, and $a_{\tau}(t)$ is the normalized synaptic
activation caused by a single pre-synaptic spike with time scale $\tau$, e.g.
$a_{\tau}(t)=e^{-t/\tau}/\tau$.

###  I.1 Continuous formulation

In the thermodynamic limit $N\to\infty$, we drop the indices in Eqs. (1), (2)
and denote $\rho(V|\eta,t)dV$ as the fraction of neurons with membrane
potentials between $V$ and $V+dV$, and parameter $\eta$ at time $t$.
Accordingly, parameter $\eta$ becomes now a continuous random variable
distributed according to a probability distribution function $g(\eta)$. The
total voltage density at time $t$ is then given by
$\int_{-\infty}^{\infty}\rho(V|\eta,t)~{}g(\eta)~{}d\eta$.

The conservation of the number of neurons leads to the following continuity
equation:

$\partial_{t}\rho+\partial_{V}\left[(V^{2}+\eta+Js+I)\rho\right]=0,$ (4)

where we have explicitly included the velocity given by (1) and (2).

##  II Results

The continuity equation (4) without temporal forcing $I(t)=0$ has a trivial
stationary solution. For each value of $\eta$, this solution has the form of a
Lorentzian function: $\rho_{0}(V|\eta)\propto(V^{2}+\eta+Js)^{-1}$.
Physically, the Lorentzian density means that firing neurons with the same
$\eta$ value will be scattered with a density inversely proportional to their
speed (1), i.e. they will accumulate at slow places and thin out at fast
places on the $V$ axis. In addition, for those $\eta$ values corresponding to
quiescent neurons, the density $\rho_{0}$ collapses at the rest state in the
form of a Dirac delta function.

Next we assume that, independently of the initial condition, solutions of (4)
generically converge to a Lorentzian-shaped function, so that all relevant
dynamics occur inside that lower dimensional space. This fact is
mathematically expressed by the following Lorentzian Ansatz (LA) for the
conditional density functions:

$\rho(V|\eta,t)=\frac{1}{\pi}\frac{x(\eta,t)}{\left[V-y(\eta,t)\right]^{2}+x(\eta,t)^{2}},$
(5)

which is a Lorentzian function with time-dependent half-width $x(\eta,t)$ and
center at $y(\eta,t)$. In the following we assume that the LA (5) completely
describes the macroscopic dynamics of the network of QIF neurons and postpone
the mathematical justification of its validity to Sec. II.5.

Figure 1: (color online). Analysis of the steady states of FREs (12). (a)
Phase diagram: In the wedge-shaped cyan-shaded region there is bistability
between a high and a low activity state. The boundary of the bistability
region is the locus of a saddle-node bifurcation which is exactly obtained in
parametric form: $(\bar{\eta},J)_{SN}=[-\pi^{2}r^{2}-3\Delta^{2}/(2\pi
r)^{2},2\pi^{2}r+\Delta^{2}/(2\pi^{2}r^{3})]$. To the right of the dashed
line, defined by $\bar{\eta}_{f}=-[J/(2\pi)]^{2}-(\pi\Delta/J)^{2}$, there is
a stable focus (shaded regions). (b) $r$-$\bar{\eta}$ and (c) $v$-$\bar{\eta}$
bifurcation diagrams for $J/\Delta^{1/2}=15$. Square symbols: Results obtained
from numerical simulations of QIF neurons (see Appendix A for details). (d)
Phase portrait of the system in the bistable region
($\bar{\eta}/\Delta=-5,J/\Delta^{1/2}=15$, triangle in panel (a)) with three
fixed points: a stable focus (with its basin of attraction shaded), a stable
node, and a saddle point.

###  II.1 Macroscopic observables: Firing rate and mean membrane potential

The half-width $x(\eta,t)$ of the LA has a particularly simple relation with
the firing rate of the neuronal population (i.e. the number of spikes per unit
time). Indeed, the firing rate for each $\eta$ value at time $t$, $r(\eta,t)$,
can be computed by noting that neurons fire at a rate given by the probability
flux at infinity:
$r(\eta,t)=\rho(V\to\infty|\eta,t)\dot{V}(V\to\infty|\eta,t)$. The limit
$V\to\infty$ on the right hand side of this equation can be evaluated within
the LA, and gives the simple identity

$x(\eta,t)=\pi r(\eta,t).$ (6)

The (total) firing rate $r(t)$ is then

$r(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}x(\eta,t)g(\eta)d\eta.$ (7)

Additionally, the quantity $y(\eta,t)$ is, for each value of $\eta$, the mean
of the membrane potential:

$y(\eta,t)=\mathrm{P.V.}\int_{-\infty}^{\infty}\rho(V|\eta,t)V\,dV.$ (8)

Here we take the Cauchy principal value, defined as
$\mathrm{P.V.}\int_{-\infty}^{\infty}f(x)dx=\lim_{R\to\infty}\int_{-R}^{R}f(x)dx$,
to avoid the otherwise ill-defined integral. The mean membrane potential is
then

$v(t)=\int_{-\infty}^{\infty}y(\eta,t)g(\eta)d\eta.$ (9)

###  II.2 Firing-rate equations

Substituting the LA (5) into the continuity equation (4), we find that, for
each value of $\eta$, variables $x$ and $y$ must obey two coupled equations
which can be written in complex form as

$\partial_{t}w(\eta,t)=i\left[\eta+Js(t)-w(\eta,t)^{2}+I(t)\right].$ (10)

where $w(\eta,t)\equiv x(\eta,t)+iy(\eta,t)$. Closing this equation requires
expressing the mean synaptic activation $s(t)$ as a function of $w(\eta,t)$.
The simplest choice is to take the limit of infinitely fast synapses, $\tau\to
0$ in (3), so that we get an equality with the firing rate: $s(t)=r(t)$. This
allows for the system of QIF neurons (1)-(3) to be exactly described by Eqs.
(10) and (7); an infinite set of integro-differential equations if $g(\eta)$
is a continuous distribution.

Equation (10) is useful for general distributions $g(\eta)$ (see Appendix B),
but a particularly sharp reduction in dimensionality is achieved if $\eta$ is
distributed according to a Lorentzian distribution of half-width $\Delta$
centered at $\bar{\eta}$:

$g(\eta)=\frac{1}{\pi}\frac{\Delta}{(\eta-\bar{\eta})^{2}+\Delta^{2}}.$ (11)

Note that this distribution accounts for the quenched variability in the
external inputs. The fact that it is Lorentzian is unrelated to the LA for the
density of membrane potentials. Assuming (11) the integrals in (7) and (9) can
be evaluated closing the integral contour in the complex $\eta$-plane and
using the residue theorem111We make an analytic continuation of $w(\eta,t)$
from real $\eta$ into complex-valued $\eta=\eta_{r}+i\eta_{i}$. This is
possible into the lower half-plane $\eta_{i}<0$, since this guarantees the
half-width $x(\eta,t)$ remains positive zero:
$\partial_{t}x(\eta,t)=-\eta_{i}>0$ at $x=0$. Therefore we closed the
integrals in (7) and (9) with an arc $|\eta|e^{i\vartheta}$ with
$|\eta|\to\infty$ and $\vartheta\in(-\pi,0)$. This contour encloses one pole
of the Lorentzian distribution (11), which written in partial fractions reads:
$g(\eta)=(2\pi
i)^{-1}[(\eta-\bar{\eta}-i\Delta)^{-1}-(\eta-\bar{\eta}+i\Delta)^{-1}]$..
Notably, the firing rate and the mean membrane potential depend only on the
value of $w$ at the pole of $g(\eta)$ in the lower half $\eta$-plane:

$\pi r(t)+iv(t)=w(\bar{\eta}-i\Delta,t).$

As a result, we only need to evaluate (10) at $\eta=\bar{\eta}-i\Delta$, and
thereby obtain a system of FREs composed of two ordinary differential
equations

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle\Delta/\pi+2rv,$ (12a)
$\displaystyle\dot{v}$ $\displaystyle=$ $\displaystyle
v^{2}+\bar{\eta}+Jr+I(t)-\pi^{2}r^{2}.$ (12b)

This nonlinear system describes the macroscopic dynamics of the population of
QIF neurons in terms of the population firing rate $r$ and mean membrane
potential $v$.

It is enlightening to compare mean-field model (12) with the equations of the
spiking neurons. Note that Eq. (12b) resembles equations (1) and (2) for the
individual QIF neuron, but _without spike resetting_. Indeed, the macroscopic
firing-rate variable $r$ enters as a negative feedback term in Eq. (12b), and
impedes the explosive growth of the mean membrane potential $v$.

This negative feedback, combined with the coupling term on the right hand side
of (12a), describe the effective interaction between the firing rate and mean
membrane potential at the network level. Therefore, the FREs (12) describe the
effect of the single-cell spike generation and reset mechanism at the network
level.

In the following we examine the dynamics described by Eq. (12), and show that
they fully reproduce the macroscopic dynamics of the network of QIF neurons,
even during episodes of strong spike synchrony.

Figure 2: (color online). The transient dynamics of an ensemble of QIF model
neurons (1)-(2) are exactly described by the FREs (12). The instantaneous
firing rate (Panels a,b) and mean membrane potential (Panels c,d) of the QIF
neurons and the FREs are depicted in black and orange, respectively. Panels
(e,f) show the raster plots of 300 randomly selected QIF neurons of the
$N=10^{4}$ in the ensemble. Left Panels (a,c,e,g): At time $t=0$, a current
$I=3$ is applied to all neurons, and set to zero again at $t=30$; stimulus
$I(t)$ shown in panel (g). Right panels (b,d,f,h): At time $t=0$ a sinusoidal
current is applied to all neurons $I(t)=I_{0}\sin(\omega t)$, with $I_{0}=3$,
$\omega=\pi/20$; stimulus $I(t)$ shown in panel (h). Parameters: $J=15$,
$\bar{\eta}=-5,\Delta=1$.

###  II.3 Analysis of the firing-rate equations

To begin with the analysis of Eq. (12), we first note that in the absence of
forcing, $I(t)=0$, the only attractors of (12) are fixed points. Figure 1(a)
shows a phase diagram of the system as a function of the mean external drive
$\bar{\eta}$ and synaptic weight $J$, both normalized by the width of the
input distribution 222The number of effective parameters can be reduced by
nondimensionalizing the system as: $\tilde{\eta}=\bar{\eta}/\Delta$,
$\tilde{J}=J/\sqrt{\Delta}$, $(\tilde{r},\tilde{v})=(r,v)/\sqrt{\Delta}$,
$\tilde{t}=\sqrt{\Delta}t$, and
$\tilde{I}(\tilde{t})=I(t/\sqrt{\Delta})/\Delta$.. There are three
qualitatively distinct regions of the phase diagram: 1 - A single stable node
corresponding to a low-activity state (white), 2 - A single stable focus
(spiral) generally corresponding to a high-activity state (gray), and 3 - A
region of bistability between low and high firing rate (cyan; see a phase
portrait of this region in Fig. 1(d)). Comparison of a sample bifurcation
diagram of the fixed points from numerical simulation of networks of QIF
neurons with that obtained from the FREs (12) shows an excellent
correspondence, see Fig. 1(b,c).

A similar phase diagram can be readily reproduced by traditional heuristic
firing-rate models, with one significant qualitative difference: the presence
of a stable focus —and hence damped oscillations. Specifically, in the gray
region of the phase diagram in Fig. 1(a), the system undergoes oscillatory
decay to the stable fixed point. This oscillatory decay occurs as well for the
high-activity state over a large extent of the region of bistability (cyan),
see e.g. Fig. 1(d).

The presence of damped oscillations at the macroscopic level reflects the
transitory synchronous firing of a fraction of the neurons in the ensemble.
This behavior is common in spiking neuron models with weak noise, and is not
captured by traditional firing rate models (see e.g. Schaffer  _et al._
(2013)).

###  II.4 Analysis of the firing-rate equations: Non-stationary inputs

To show that the FREs (12) fully describe the macroscopic response of the
population of QIF neurons to time-varying stimuli (up to finite-size effects),
we consider two types of stimulus $I(t)$: 1 - a step function and 2 - a
sinusoidal forcing. In both cases we simulate the full system of QIF neurons
and the FREs (12).

Figure 2 shows the system’s response to the two different inputs. In both
cases the system is initially ($t<0$) in a bistable regime and set in the low
activity state, with parameters corresponding to those of Fig. 1(d). Left
panels of Fig. 2 show the response of the system to a step current, applied at
$t=0$. The applied current is such that the system abandons the bistable
region —see Fig. 1(a)— and approaches the high activity state, which is a
stable focus. This is clearly reflected in the time series $r(t),~{}v(t)$,
where the rate equations (12) exactly predict the damped oscillations
exhibited by the mean-field of the QIF neurons. The raster plot in panel (e)
shows the presence of the oscillations, which is due to the transitory
synchronous firing of a large fraction of neurons in the population. Finally,
at $t=30$, the current is removed and the system converges —again, showing
damped oscillations— to the new location of the (focus) fixed point, which
clearly coexists with the stable node where it was originally placed ($t<0$).

The right hand panels of Fig. 2 show the response of the model to a periodic
current, which drives the system from one side of the bistable region to the
other. As a result, we observe periodic bursting behavior when the system
visits the stable focus region of the phase diagram Fig. 1(a).

To further illustrate the potential of the FREs (12) to predict and
investigate complex dynamics in ensembles of spiking neurons, we present here
a simple situation where the system of QIF neurons exhibits macroscopic chaos.
This is observed by increasing the frequency $\omega$ of the sinusoidal
driving, so that the system cannot trivially follow the stable fixed point at
each cycle of the applied current.

Figure 3: (color online). Firing rate model predicts the existence of
macroscopic chaos in a ensemble QIF neurons (1) with an injected periodic
current $I(t)=I_{0}\sin(\omega t)$. (a) Phase diagram with the regions where
chaos is found, obtained using the rate model (12). In the black-shaded
region, there is only one chaotic attractor, whereas in the cyan region the
chaotic attractor coexists with a periodic attractor. Dotted lines correspond
to the bistability boundary of Fig. 1(a), depicted to facilitate comparison.
(b) Chaotic trajectory obtained simulating the FREs (12); the Lyapunov
exponent is $\lambda=0.183\ldots$. (c) Chaotic trajectory obtained from the
QIF neurons —parameters corresponding to the (yellow) triangle symbol in panel
(a). (d) Time series for the rate model (orange) and QIF model (black). (e)
Raster plot corresponding to 300 randomly selected neurons. (f) Same raster
plot as in (e), ordering neurons according to their intrinsic current
$\eta_{k}$. Parameters: $I_{0}=3$, $\omega=\pi$; and $\bar{\eta}=-2.5$,
$J=10.5$ (triangle symbol in (a)).

Figure 3(a) shows a phase diagram obtained using Eq. (12). The shaded regions
indicate parameter values where the rate model has either a single chaotic
attractor (in black) or a chaotic attractor coexisting with a periodic orbit
(in cyan). A trajectory on the chaotic attractor is depicted in Fig. 3(b), and
the clearly irregular time series of the firing rate is shown in orange in
Fig. 3(d). Using the same parameters we performed numerical simulations of the
QIF neurons (1)-(2), finding a similar attractor and irregular dynamics as in
the rate model, see Figs. 3(c) and (d). To obtain the time series shown in
Figs. 3(d), we ran the QIF neurons numerically and, after a long transient, at
time $t=0$, the rate model (12) was initiated with the values of $r$ and $v$
obtained from the population of QIF neurons. The time series of the two
systems, which are initially close, rapidly diverge reflecting the chaotic
nature of the system.

Finally, to illustrate this chaotic state at the microscopic level, Fig. 3(e)
shows the raster plot for 300 randomly chosen neurons, corresponding to the
time series in Fig. 3(d). The irregular firing of neurons in Fig. 3(e) has
some underlying structure that may be visualized ordering the same set of
neurons according to their intrinsic currents as $\eta_{k}<\eta_{k+1}$, with
$k\in[1,300]$, see Fig. 3(f). Clearly, the maxima of the firing rate coincide
with the synchronous firings of clusters of neurons with similar $\eta$
values. The size of these clusters is highly irregular in time, in
concomitance with the chaotic behavior.

###  II.5 Validity of the Lorentzian ansatz

Thus far, we have shown that the LA (5) solves the continuity equation (4),
and confirmed that these solutions agree with the numerical simulations of the
original system of QIF neurons (1)-(2). Here we further clarify why the LA
holds for ensembles of QIF neurons.

Transforming the voltage of the QIF neuron into a phase via
$V_{j}=\tan(\theta_{j}/2)$, the system (1)-(2) becomes an ensemble of ‘theta
neurons’ Ermentrout and Kopell (1986):

$\dot{\theta}_{j}=(1-\cos\theta_{j})+(1+\cos\theta_{j})[\eta_{j}+Js(t)+I(t)].$
(13)

In the new phase variable $\theta\in[0,2\pi)$, the LA (5) becomes:

$\tilde{\rho}(\theta|\eta,t)=\frac{1}{2\pi}\mbox{Re}\left[\frac{1+\alpha(\eta,t)e^{i\theta}}{1-\alpha(\eta,t)e^{i\theta}}\right],$
(14)

where the function $\alpha(\eta,t)$ is related with $w(\eta,t)$ as

$\alpha(\eta,t)=\frac{1-w(\eta,t)}{1+w(\eta,t)}.$ (15)

Equations (5) and (14) are two representations of the so-called Poisson kernel
on the half-plane and on the unit disk, respectively. These representations
are related via equation (15), that establishes a conformal mapping from the
half-plane $\mbox{Re}(w)\geq 0$ onto the unit disk $|\alpha|\leq 1$. In the
next section we show that variables $r$ and $v$ can be related, via the same
conformal map, with the Kuramoto order parameter, which is a macroscopic
measure of phase coherence Kuramoto (1984); Pikovsky  _et al._ (2001).

The key observation supporting the applicability of the LA is the fact that
equation (14) turns out to be the ansatz discovered in 2008 by Ott and
Antonsen Ott and Antonsen (2008). According to the Ott-Antonsen (OA) theory,
in the thermodynamic limit the dynamics of the class of systems

$\partial_{t}\theta(\eta,t)=\Omega(\eta,t)+\mathrm{Im}\left[H(\eta,t)e^{-i\theta}\right],$
(16)

generally converges to the OA manifold (14). In our case, for equation (13),
we have $\Omega(\eta,t)=1+\eta+Js+I$ and $H(\eta,t)=i(-1+\eta+Js+I)$. Thus
far, the convergence of (16) to the OA manifold has been proven only for
$H(\eta,t)=H(t)$. This includes the well-known Kuramoto and Winfree models Ott
and Antonsen (2009); *OHA11; Pazó and Montbrió (2014), but not system (13).
However, there are theoretical arguments Pikovsky and Rosenblum (2011) that
strongly suggest that the OA manifold is also attracting for $\eta$-dependent
$H$, as numerically confirmed by a number of recent papers using theta neurons
Luke  _et al._ (2013); So  _et al._ (2014); Laing (2014) and other phase-
oscillator models Montbrió and Pazó (2011a); Pazó and Montbrió (2011);
Montbrió and Pazó (2011b); Iatsenko  _et al._ (2013, 2014).

###  II.6 Firing rate and Kuramoto order parameter

There exists a mapping between the macroscopic variables $r$ and $v$ and the
so-called Kuramoto order parameter $Z$. Equation (15) relates, for each value
of $\eta$, the firing rate and the mean membrane potential, both contained in
$w$, with the uniformity of the phase density, measured by $\alpha$ —note that
$\alpha=0$ in (14) yields a perfectly uniform density. The Kuramoto order
parameter is obtained by integrating $\alpha$ over the whole population as

$\displaystyle Z(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g(\eta)\int_{0}^{2\pi}\tilde{\rho}(\theta|\eta,t)e^{i\theta}d\theta\,d\eta$
(17) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g(\eta)\alpha^{*}(\eta,t)\,d\eta,$

where we assumed that the density $\tilde{\rho}$ is the OA ansatz Eq. (14).
The quantity $Z$ can be seen as the center of mass of a population of phases
distributed across the unit circle $e^{i\theta}$.

For a Lorentzian distribution of currents, $g(\eta)$, we may derive the exact
formula (see Appendix B)

$Z=\frac{1-W^{*}}{1+W^{*}}.$ (18)

relating the Kuramoto order parameter $Z$ with the firing-rate quantity
$W\equiv\pi r+iv$. Figure 4 illustrates this conformal mapping of the right
half-plane ($r\geq 0$) onto the unit disk ($|Z|\leq 1$).

Figure 4: (color online). The conformal map (18) transforms the right half-
plane onto the unit disk. This transformation and its inverse define a one-to-
one mapping between the Kuramoto order parameter $Z$ and the macroscopic
quantities of the population of QIF neurons: $r$: firing rate, and $v$: mean
membrane potential. Note that in the limit $Z=e^{i\Psi}$ (full coherence,
$\tilde{\rho}(\theta)=\delta(\theta-\Psi)$) we recover
$v=i(1-e^{i\Psi})/(1+e^{i\Psi})=\tan(\Psi/2)$, as the original mapping
$V=\tan(\theta/2)$ between QIF- and theta-neurons dictates.

##  III Conclusions

We have presented a method for deriving firing rate equations for a network of
heterogeneous QIF neurons, which is exact in the thermodynamic limit. To our
knowledge, the resulting system of ordinary differential equations (12)
represents the first exact FREs of a network of spiking neurons.

We would like to emphasize that the derivation of the FREs does not rely on
assumptions of weak coupling, separation of time-scales, averaging, or any
other approximation. Rather, the assumptions underlying the validity of
equation (12) are that $i)$ the QIF neurons be all-to-all connected, $ii)$
heterogeneity be quenched, and $iii)$ inputs be distributed according to a
Lorentzian distribution.

The last two assumptions may seem particularly restrictive. Concerning
$(iii)$, it must be stressed that for _arbitrary_ distributions of quenched
heterogeneity the LA (5) remains generally valid, and therefore equation (10)
can be used to discern the stability of the network states (see Appendix C).
Furthermore, in Appendix C we also show numerical simulations using uniform-
and Gaussian-distributed inputs, that reveal very similar macroscopic dynamics
in response to time varying inputs. Even relaxing assumption $(ii)$, numerical
simulations with identical QIF neurons driven by external independent Gaussian
noise sources show qualitative agreement with the FREs (12) (see Appendix D).
In sum, the choice of a quenched Lorentzian distribution $g(\eta)$ is thus a
mere mathematical convenience, whereas the insights gained from the resulting
firing-rate description are valid more generally.

Our derivation represents a sharp departure from, and a major advance compared
to previous studies in several regards. Firstly, in the past it has only been
possible to calculate the approximate firing rate of networks of spiking
neurons for stationary states, or for weak deviations of such states Gerstner
(1995); Fusi and Mattia (1998); Amit and Brunel (1997); Brunel (2000); Shriki
_et al._ (2003); Ostojic and Brunel (2011); Schaffer  _et al._ (2013). The
equations which result from these derivations are difficult to solve, and
typically require special numerical methods. In contrast, FREs (12) can be
easily analyzed and simulated, and exactly reproduce the behavior of the
spiking network far from any fixed point and for arbitrary external currents.
Secondly, firing-rate descriptions traditionally assume that the activity of a
population of neurons is equivalent to a set of uncorrelated stochastic
processes with a given rate, see e.g. Dayan and Abbott (2001); Gerstner and
Kistler (2002). However, in simulations of spiking networks it is well known
that the response of the network to non-stationary inputs generally involves
some degree of spike synchrony. Recent theoretical work has sought to improve
on classical rate models, which act as a low-pass filter of inputs, by fitting
them with linear filters extracted from the corresponding Fokker-Planck
equation for the network Ostojic and Brunel (2011); Schaffer  _et al._ (2013).
These filters tend to generate damped oscillations when the external noise
level is not too high, reflecting the presence of spike synchrony Shriki  _et
al._ (2003); Schaffer  _et al._ (2013). The FREs (12) also capture this
phenomenon, and reveal that the underlying mechanism is an interplay between
firing rate and subthreshold voltage. Furthermore, because these equations are
exact, we can explore the full nonlinear response of the network, such as the
generation of chaotic states of synchronous bursting as shown in Fig. 3.

Recently it has been shown that, not only Kuramoto-like models, but a much
wider class of phase-models, evolve in the OA manifold defined by (14).
Specifically, networks of pulse-coupled oscillators Pazó and Montbrió (2014)
and theta-neurons Luke  _et al._ (2013); So  _et al._ (2014); Laing (2014)
allow for an exact, low-dimensional description in terms of the Kuramoto order
parameter (17). Although these later works use a finite-width phase coupling
function that differs from the synaptic coupling (3), the obtained low-
dimensional description in terms of the Kuramoto order parameter is analogous
to ours, but in a different space 333However, note that the phase diagrams in
Luke et al. Luke  _et al._ (2013) do not fully agree with ours, and show
additional instabilities. These bifurcations are associated with the finite
width of the pulses, as it accurs in the Winfree model Pazó and Montbrió
(2014).. Indeed, we showed that the OA ansatz (14) is related, via the
nonlinear transformation of variables (15), to the LA ansatz (5). Remarkably,
this transformation establishes an exact correspondence between the Kuramoto
order parameter and a novel, biophysically maningful macroscopic observable
which describes the firing rate and mean membrane potential of the neuronal
network. Interestingly, the low-dimensional description in terms of firing
rates seems to be a more natural description for networks of spiking neurons,
compared to that in terms of the Kuramoto order parameter. The firing rate
equations (12) take a surprisingly simple form in the LA manifold, what makes
them a valuable tool to explore and and understand the mechanisms governing
the macroscopic dynamics of neuronal networks.

Finally, since the OA ansatz is the asymptotic solution for systems of the
form in equation (16), applying the change of variables $V=\tan(\theta/2)$
automatically implies that the LA should hold for populations governed by:

$\dot{V}_{j}=A(\eta_{j},t)+B(\eta_{j},t)V_{j}+C(\eta_{j},t)V_{j}^{2},$ (19)

where $A$, $B$, and $C$ are related with $\Omega$ and $H$ as
$A=[\Omega+\mathrm{Im}(H)]/2$, $B=-\mathrm{Re}(H)$,
$C=[\Omega-\mathrm{Im}(H)]/2$. Notably, equation (19) defines a wide family of
ensembles of QIF neurons.

Therefore, the LA is actually valid for more general networks beyond the one
investigated here —$\eta_{j}$ in Eq. (19) may also be a vector containing
several forms of disorder. As particularly relevant cases, in Appendix E we
provide the FREs governing the dynamics of an excitatory network with both,
heterogeneous inputs $\eta$ _and_ synaptic weights $J$, as well as pair of
interacting excitatory and inhibitory populations of QIF neurons. In addition,
according to Eq. (19), the LA is also valid if synapses are modeled as
conductances, in which case the reversal potentials may be distributed as
well. Moreover, the role of gap-junctions or synaptic kinetics may be
considered within the same framework.

## References

  * Averbeck  _et al._ (2006) B. B. Averbeck, P. E. Latham, and A. Pouget, “Neural correlations, population coding and computation,” Nature Reviews Neuroscience 7, 358–366 (2006).
  * Hodgkin and Huxley (1952) A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” J Physiol. 117, 500–544 (1952).
  * Brette  _et al._ (2007) Romain Brette, Michelle Rudolph, Ted Carnevale, Michael Hines, David Beeman, James M Bower, Markus Diesmann, Abigail Morrison, Philip H Goodman, Frederick C Harris Jr, _et al._ , “Simulation of networks of spiking neurons: a review of tools and strategies,” Journal of computational neuroscience 23, 349–398 (2007).
  * Izhikevich and Edelman (2008) E. M. Izhikevich and G. M. Edelman, “Large-scale model of mammalian thalamocortical systems,” Proc. Nat. Acad. Sci. 105, 3593–3598 (2008).
  * Izhikevich (2007) E. M. Izhikevich,  _Dynamical Systems in Neuroscience._ (The MIT Press, Cambridge, Massachusetts, 2007).
  * Wilson and Cowan (1972) Hugh R Wilson and Jack D Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons,” Biophys. J. 12, 1–24 (1972).
  * Amari (1974) Shun-ichi Amari, “A method of statistical neurodynamics,” Kybernetik 14, 201–215 (1974).
  * Freeman (1975) Walter J Freeman,  _Mass action in the nervous system_ (Academic Press, New York, 1975).
  * Gerstner (1995) W. Gerstner, “Time structure of the activity in neural network models,” Phys. Rev. E 51, 738–758 (1995).
  * Fusi and Mattia (1998) S. Fusi and M. Mattia, “Collective behavior of networks with linear (vlsi) integrate-and-fire neurons,” Neural Comput. 11, 633–652 (1998).
  * Amit and Brunel (1997) Daniel J Amit and Nicolas Brunel, “Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex.” Cerebral Cortex 7, 237–252 (1997).
  * Brunel (2000) N. Brunel, “Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons,” Journal of computational neuroscience 8, 183–208 (2000).
  * Dayan and Abbott (2001) P. Dayan and L. F. Abbott,  _Theoretical neuroscience_ (Cambridge, MA: MIT Press, 2001).
  * Gerstner and Kistler (2002) W. Gerstner and W. M. Kistler,  _Spiking neuron models: Single neurons, populations, plasticity_ (Cambridge university press, 2002).
  * Shriki  _et al._ (2003) Oren Shriki, David Hansel, and Haim Sompolinsky, “Rate models for conductance-based cortical neuronal networks,” Neural Comput. 15, 1809–1841 (2003).
  * Coombes (2005) Stephen Coombes, “Waves, bumps, and patterns in neural field theories,” Biol. Cybern. 93, 91–108 (2005).
  * Roxin  _et al._ (2005) A. Roxin, N. Brunel, and D. Hansel, “Role of delays in shaping spatiotemporal dynamics of neuronal activity in large networks,” Phys. Rev. Lett. 94, 238103 (2005).
  * Deco  _et al._ (2008) G. Deco, V. K. Jirsa, P. A. Robinson, M. Breakspear, and K. Friston, “The dynamic brain: from spiking neurons to neural masses and cortical fields,” PLoS Computational biology 4, e1000092 (2008).
  * Ermentrout and Terman (2010) G. B. Ermentrout and D. H. Terman,  _Mathematical foundations of neuroscience_ , Vol. 64 (Springer, 2010).
  * Ostojic and Brunel (2011) S. Ostojic and N. Brunel, “From spiking neuron models to linear-nonlinear models,” PLoS Comput. Biol. 7, e1001056 (2011).
  * Roxin and Montbrió (2011) A. Roxin and E. Montbrió, “How effective delays shape oscillatory dynamics in neuronal networks,” Physica D 240, 323–345 (2011).
  * Schaffer  _et al._ (2013) Evan S Schaffer, Srdjan Ostojic, and LF Abbott, “A complex-valued firing-rate model that approximates the dynamics of spiking networks,” PLoS Comput. Biol. 9, e1003301 (2013).
  * Hopfield (1984) J. J. Hopfield, “Neurons with graded response have collective computational properties like those of two-state neurons,” Proceedings of the national academy of sciences 81, 3088–3092 (1984).
  * Mongillo  _et al._ (2008) G. Mongillo, O. Barak, and M. Tsodyks, “Synaptic theory of working memory,” Science 319, 1543–1546 (2008).
  * Ben-Yishai  _et al._ (1995) R Ben-Yishai, R Lev Bar-Or, and H Sompolinsky, “Theory of orientation tuning in visual cortex.” Proc. Nat. Acad. Sci. 92, 3844–3848 (1995).
  * Hansel and Sompolinsky (1998) David Hansel and Haim Sompolinsky, “Modeling feature selectivity in local cortical circuits,” in  _Methods in Neuronal Modelling: From Ions to Networks_ , edited by C. Koch and I. Segev (MIT Press, Cambridge, 1998) pp. 499–567.
  * Moreno-Bote  _et al._ (2007) Rubén Moreno-Bote, John Rinzel, and Nava Rubin, “Noise-induced alternations in an attractor network model of perceptual bistability,” J. Neurophysiol. 98, 1125–1139 (2007).
  * Zhang (1996) K. Zhang, “Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory,” J. Neurosci. 16, 2112–2126 (1996).
  * Wong and Wang (2006) K.-F. Wong and X.-J. Wang, “A recurrent network mechanism of time integration in perceptual decisions,” J. Neurosci. 26, 1314–1328 (2006).
  * Abeles (1991) Moshe Abeles,  _Corticonics: Neural circuits of the cerebral cortex_ (Cambridge University Press, 1991).
  * Varela  _et al._ (2001) F. Varela, J.-P. Lachaux, E. Rodriguez, and J. Martinerie, “The brainweb: phase synchronization and large-scale integration,” Nature reviews neuroscience 2, 229–239 (2001).
  * Engel and Singer (2001) Andreas K Engel and Wolf Singer, “Temporal binding and the neural correlates of sensory awareness,” Trends in cognitive sciences 5, 16–25 (2001).
  * Ward (2003) L. M. Ward, “Synchronous neural oscillations and cognitive processes,” Trends in cognitive sciences 7, 553–559 (2003).
  * Fries (2005) P. Fries, “A mechanism for cognitive dynamics: neuronal communication through neuronal coherence,” Trends in cognitive sciences 9, 474–480 (2005).
  * Buzsaki (2006) Gyorgy Buzsaki,  _Rhythms of the Brain_ (Oxford University Press, 2006).
  * Fell and Axmacher (2011) J. Fell and N. Axmacher, “The role of phase synchronization in memory processes,” Nature Reviews Neuroscience 12, 105–118 (2011).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” Chaos 18, 037113 (2008).
  * Hodgkin (1948) AL Hodgkin, “The local electric changes associated with repetitive action in a non-medullated axon,” The Journal of physiology 107, 165–181 (1948).
  * Rinzel and Ermentrout (1998) J. Rinzel and B. Ermentrout, “Analysis of neural excitability and oscillations,” in  _Methods in Neuronal Modelling: From Ions to Networks_ , edited by C. Koch and I. Segev (MIT Press, Cambridge, 1998) pp. 251–291.
  * Ermentrout and Kopell (1986) B. Ermentrout and N. Kopell, “Parabolic bursting in an excitable system coupled with a slow oscillation,” SIAM J. Appl. Math. 46, 233–253 (1986).
  * Latham  _et al._ (2000) PE Latham, BJ Richmond, PG Nelson, and S Nirenberg, “Intrinsic dynamics in neuronal networks. i. theory,” Journal of Neurophysiology 83, 808–827 (2000).
  * Buice and Chow (2013) M. A. Buice and C. C. Chow, “Dynamic finite size effects in spiking neural networks,” PLoS Comput. Biol. 9, e1002872 (2013).
  * Note (1) We make an analytic continuation of $w(\eta,t)$ from real $\eta$ into complex-valued $\eta=\eta_{r}+i\eta_{i}$. This is possible into the lower half-plane $\eta_{i}<0$, since this guarantees the half-width $x(\eta,t)$ remains positive zero: $\partial_{t}x(\eta,t)=-\eta_{i}>0$ at $x=0$. Therefore we closed the integrals in (7) and (9) with an arc $|\eta|e^{i\vartheta}$ with $|\eta|\to\infty$ and $\vartheta\in(-\pi,0)$. This contour encloses one pole of the Lorentzian distribution (11), which written in partial fractions reads: $g(\eta)=(2\pi i)^{-1}[(\eta-\mathaccentV{bar}016\eta-i\Delta)^{-1}-(\eta-\mathaccentV{bar}016\eta+i\Delta)^{-1}]$.
  * Note (2) The number of effective parameters can be reduced by nondimensionalizing the system as: $\mathaccentV{tilde}07E\eta=\mathaccentV{bar}016\eta/\Delta$, $\mathaccentV{tilde}07EJ=J/\sqrt{\Delta}$, $(\mathaccentV{tilde}07Er,\mathaccentV{tilde}07Ev)=(r,v)/\sqrt{\Delta}$, $\mathaccentV{tilde}07Et=\sqrt{\Delta}t$, and $\mathaccentV{tilde}07EI(\mathaccentV{tilde}07Et)=I(t/\sqrt{\Delta})/\Delta$.
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer-Verlag, Berlin, 1984).
  * Pikovsky  _et al._ (2001) A. S. Pikovsky, M. G. Rosenblum, and J. Kurths,  _Synchronization, a Universal Concept in Nonlinear Sciences_ (Cambridge University Press, Cambridge, 2001).
  * Ott and Antonsen (2009) E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator systems,” Chaos 19, 023117 (2009).
  * Ott  _et al._ (2011) E. Ott, B. R. Hunt, and T. M. Antonsen, “Comment on “long time evolution of phase oscillators systems”,” Chaos 21, 025112 (2011).
  * Pazó and Montbrió (2014) D. Pazó and E. Montbrió, “Low-dimensional dynamics of populations of pulse-coupled oscillators,” Phys. Rev. X 4, 011009 (2014).
  * Pikovsky and Rosenblum (2011) A. Pikovsky and M. Rosenblum, “Dynamics of heterogeneous oscillator ensembles in terms of collective variables,” Physica D 240, 872 – 881 (2011).
  * Luke  _et al._ (2013) T. B. Luke, E. Barreto, and P. So, “Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons,” Neural Comput. 25, 3207–3234 (2013).
  * So  _et al._ (2014) P. So, T. B. Luke, and E. Barreto, “Networks of theta neurons with time-varying excitability: Macroscopic chaos, multistability, and final-state uncertainty,” Physica D 267, 16–26 (2014).
  * Laing (2014) C. R. Laing, “Derivation of a neural field model from a network of theta neurons,” Phys. Rev. E 90, 010901 (2014).
  * Montbrió and Pazó (2011a) E. Montbrió and D. Pazó, “Shear diversity prevents collective synchronization,” Phys. Rev. Lett. 106, 254101 (2011a).
  * Pazó and Montbrió (2011) D. Pazó and E. Montbrió, “The Kuramoto model with distributed shear,” EPL (Europhys. Lett.) 95, 60007 (2011).
  * Montbrió and Pazó (2011b) E. Montbrió and D. Pazó, “Collective synchronization in the presence of reactive coupling and shear diversity,” Phys. Rev. E 84, 046206 (2011b).
  * Iatsenko  _et al._ (2013) D. Iatsenko, S. Petkoski, P. V. E. McClintock, and A. Stefanovska, “Stationary and traveling wave states of the Kuramoto model with an arbitrary distribution of frequencies and coupling strengths,” Phys. Rev. Lett. 110, 064101 (2013).
  * Iatsenko  _et al._ (2014) D. Iatsenko, P. V. E. McClintock, and A. Stefanovska, “Oscillator glass in the generalized Kuramoto model: synchronous disorder and two-step relaxation,” Nat. Commun. 5, 4188 (2014).
  * Note (3) However, note that the phase diagrams in Luke et al. Luke  _et al._ (2013) do not fully agree with ours, and show additional instabilities. These bifurcations are associated with the finite width of the pulses, as it accurs in the Winfree model Pazó and Montbrió (2014).
  * Roberts and Quispel (1992) J.A.G. Roberts and G.R.W. Quispel, “Chaos and time-reversal symmetry. Order and chaos in reversible dynamical systems,” Phys. Rep. 216, 63 – 177 (1992).
  * Strogatz and Mirollo (1991) S. H. Strogatz and R. E. Mirollo, “Stability of incoherence in a population of coupled oscillators,” J. Stat. Phys. 63, 613–635 (1991).
  * Strogatz  _et al._ (1992) S. H. Strogatz, R. E. Mirollo, and P. C. Matthews, “Coupled nonlinear oscillators below the synchronization threshold: Relaxation by generalized Landau damping,” Phys. Rev. Lett. 68, 2730–2733 (1992).
  * Chiba and Nishikawa (2011) H. Chiba and I. Nishikawa, “Center manifold reduction for large populations of globally coupled phase oscillators,” Chaos 21, 043103 (2011).
  * Daido (1996) H. Daido, “Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions: bifurcation of the order function,” Physica D 91, 24–66 (1996).
  * Note (4) Equation (10) becomes then a reversible system Roberts and Quispel (1992) of ordinary differential equations with the invariance under time reversal and $w\to w^{*}$, e.g. Eq. (12) for $\Delta=0$ is reversible, with invariance under $t\to-t$, $v\to-v$. It is advisable to approximate Eq. (10) by a set of ODEs (one for each $\eta$ value) with additional small terms breaking time-reverisibility. This is achieved including a small term in each ODE, like $\Delta$ in Eq. (12a), so that $g(\eta)$ is actually approximated by a sum of very sharp Lorentzian functions.
  * Watanabe and Strogatz (1994) S. Watanabe and S. H. Strogatz, “Constant of motion for superconducting Josephson arrays,” Physica D 74, 197–253 (1994).
  * Pikovsky and Rosenblum (2008) A. Pikovsky and M. Rosenblum, “Partially integrable dynamics of hierarchical populations of coupled oscillators,” Phys. Rev. Lett. 101, 264103 (2008).
  * Note (5) A subtle question is that, as it occurs in the Kuramoto model, there exists a continuous spectrum of eigenvalues at $\lambda=-2w_{0}(\eta)i$ and the complex conjugate. Indeed, for $\eta$ values such that neurons are actively firing ($y_{0}(\eta)=0$), the eigenvalues lie exactly on the imaginary axis. Apparently, this seems to forbid the exponential stability obtained from Eq. (12) for the Lorentzian distribution. The same mathematical conundrum was encountered and solved in the Kuramoto model Strogatz and Mirollo (1991); Strogatz  _et al._ (1992); Chiba and Nishikawa (2011), what suggests similar mechanisms are operating here.
  * Note (6) Note that considering stronger forms of heterogeneity for the synaptic weights may generally not lead to Eq. (19).

## Appendix A: Numerical Simulations

To numerically simulate the population of QIF neurons (1) we used the Euler
method, with time step $dt=10^{-4}$. The population had $N=10^{4}$ neurons,
and the Lorentzian distribution (11) was deterministically generated using:
$\eta_{j}=\bar{\eta}+\Delta\arctan[\pi/2(2j-N-1)/(N+1)]$, where $j=1,\dots,N$,
and $\Delta=1$.

The time it takes for the membrane potential of a QIF neuron (with $I_{j}>0$)
to reach infinity from a given positive value of the membrane potential ̵̇is
$\arctan(\sqrt{I_{j}}/V_{p})/\sqrt{I_{j}}$. For $\sqrt{I_{j}}\ll V_{p}$, this
expression can be approximated as:

$\frac{\arctan\left(\sqrt{I_{j}}/V_{p}\right)}{\sqrt{I_{j}}}\approx\frac{1}{V_{p}}.$

In simulations we considered $V_{p}=-V_{r}=100$, and used the previous
approximation. Thus, the time for the neurons to reach infinity from $V_{p}$
is $1/V_{p}\approx 10^{-2}$. Additionally, the time taken from minus infinity
to $V_{r}$ is $−1/V_{r}\approx 10^{-2}$. Numerically, once a neuron’s membrane
potential $V_{j}$ satisfies $V_{j}\geq V_{p}$ the neuron is reset to $-V_{j}$
and held there for a refractory time $2/V_{j}$. The neurons produce a spike
when $V_{j}$ reaches infinity, i.e. a time interval $1/V_{j}$ after crossing
$V_{p}$. The exact time of the spike can not be evaluated exactly in numerical
simulations, since $dt$ is finite. However, simulations agree very well with
the theory provided that $1/V_{p}\gg dt$.

To evaluate the mean membrane potential $v$, the population average was
computed discarding those neurons in the refractory period. The choice
$V_{p}=-V_{r}$ is not needed, but in this way the observed mean membrane
potential $v$ agrees with the theory, consistent with the definition in (8).
The mean synaptic activation (3) was evaluated using the Heaviside step
function $a_{\tau}(t)=\Theta(\tau-t)/\tau$ with $\tau=10^{-3}$. The
instantaneous firing rates were obtained binning time and counting spikes
within a sliding time window of size $\delta t=2\times 10^{-2}$ (Figs. 2, C2,
C3, D1) and $\delta t=4\times 10^{-2}$ (Fig. 3).

## Appendix B: Proof of equation (18)

The inverse of Eq. (18) reads:

$W=\frac{1-Z^{*}}{1+Z^{*}}$ (B1)

Recalling Eqs. (7) and (9), we write the macroscopic field $W$ as

$W(t)=\int_{-\infty}^{\infty}w(\eta,t)g(\eta)d\eta.$

where $w\equiv x+iy$. Inserting the conformal mapping
$w=(1-\alpha)/(1+\alpha)$ —the inverse of Eq. (15)– we get

$W(t)=\int_{-\infty}^{\infty}\left[\frac{1-\alpha(\eta,t)}{1+\alpha(\eta,t)}\right]g(\eta)d\eta.$

Using the geometric series formula, and grouping the powers of $\alpha$, we
may evaluate the integrals obtaining

$W(t)=1-2Z^{*}(t)+2Z_{2}^{*}(t)-2Z_{3}^{*}(t)+\cdots,$ (B2)

where the $Z_{m}$ are the generalized order parameters Daido (1996):

$\displaystyle Z_{m}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g(\eta)\int_{0}^{2\pi}\tilde{\rho}(\theta|\eta,t)e^{im\theta}d\theta\,d\eta$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g(\eta)[\alpha^{*}(\eta,t)]^{m}\,d\eta.$

Since, for Lorentzian $g(\eta)$,
$Z_{m}(t)=[\alpha^{*}(\eta=\bar{\eta}-i\Delta,t)]^{m}=Z(t)^{m}$, we can revert
the power series in (B2) and obtain the result in (B1).

## Appendix C: Results for arbitrary distributions of currents

In this Appendix we compare the results in the main text, obtained using a
Lorentzian distribution of currents $g(\eta)$, with other distributions. The
results are qualitatively similar as evidenced by Figs. C1, C2 and C3 (cf.
Figs. 1, 2, and 3 in the main text).

Note that if a particular distribution $g(\eta)$ has $2n$ poles (all of them
off the real axis) one can readily obtain the FREs consisting of $n$ complex-
valued ordinary differential equations. Even if $g(\eta)$ does not fulfill
this condition, as in the case of Gaussian or uniform distributions, the
Lorentzian Ansatz still holds (see below).

Moreover, it is possible to efficiently simulate the dynamics of a population
by integrating Eq. (10) for a sample of $\eta$ values that approximate a
particular $g(\eta)~{}$444Equation (10) becomes then a reversible system
Roberts and Quispel (1992) of ordinary differential equations with the
invariance under time reversal and $w\to w^{*}$, e.g. Eq. (12) for $\Delta=0$
is reversible, with invariance under $t\to-t$, $v\to-v$. It is advisable to
approximate Eq. (10) by a set of ODEs (one for each $\eta$ value) with
additional small terms breaking time-reverisibility. This is achieved
including a small term in each ODE, like $\Delta$ in Eq. (12a), so that
$g(\eta)$ is actually approximated by a sum of very sharp Lorentzian
functions.. However, if $g(\eta)$ is a discrete distribution (or has a
discrete part), —as in the case of the Dirac delta function
$g(\eta)=\delta(\eta-\bar{\eta})$—, the so-called Watanabe-Strogatz theory
applies Watanabe and Strogatz (1994); Pikovsky and Rosenblum (2011, 2008), and
the ‘Lorentzian manifold’ is no longer attracting since the dynamics is
extremely degenerate, akin to a Hamiltonian system.

### Steady states

Considering $I(t)=0$ in Eqs. (1) and (2), it is possible to find the solutions
with steady firing rate $r(t)=r_{0}$ for arbitrary distributions of currents.

As we pointed out in the main text, the stationary solution of the continuity
equation (4) for those neurons that are intrinsically active ($\eta>-Jr_{0}$)
must be is inversely proportional to their speed, that is

$\rho_{0}(V|\eta)=\frac{c(\eta)}{V^{2}+\eta+Jr_{0}},$

where the normalizing constant is $c(\eta)=\sqrt{\eta+Jr_{0}}/\pi$. The firing
rate for each value of $\eta$ is then

$r_{0}(\eta)=\rho_{0}(V\to\infty|\eta)\dot{V}(V\to\infty|\eta,t)=c(\eta),$

and integration over all the firing neurons, gives the self-consistent
condition for $r_{0}$:

$r_{0}=\int_{-Jr_{0}}^{\infty}\,c(\eta)g(\eta)\,d\eta$ (C1)

which is valid for any distribution $g(\eta)$ (an analogous expression was
obtained in Eq. (45) of Buice and Chow (2013)). For Lorentzian $g(\eta)$, the
integral in (C1) can be evaluated, and solving the resulting equation for
$r_{0}$ one obtains the steady states in agreement with the result of setting
$\dot{r}=\dot{v}=0$ in the FREs (12).

### Linear stability analysis of steady states

Assuming that the LA captures the actual dynamics of the population of QIF
neurons, one can also investigate the stability of steady states for arbitrary
distributions $g(\eta)$. Indeed, the evolution of infinitesimal perturbations
from a steady state

$w_{0}(\eta)\equiv
x_{0}(\eta)+iy_{0}(\eta)=\begin{cases}\sqrt{\eta+Jr_{0}}&\text{if
}\eta>-Jr_{0}\\\ -i\sqrt{-\eta-Jr_{0}}&\text{if }\eta\leq-Jr_{0}\end{cases}$

is determined by the linearization of Eq. (10):

$\partial_{t}\delta w(\eta,t)=i(J\delta r-2w_{0}\delta w)$ (C2)

where $\delta r(t)=(2\pi)^{-1}\int_{-\infty}^{\infty}[\delta w(\eta,t)+\delta
w^{*}(\eta,t)]g(\eta)d\eta$.

To find the discrete spectrum of eigenvalues 555A subtle question is that, as
it occurs in the Kuramoto model, there exists a continuous spectrum of
eigenvalues at $\lambda=-2w_{0}(\eta)i$ and the complex conjugate. Indeed, for
$\eta$ values such that neurons are actively firing ($y_{0}(\eta)=0$), the
eigenvalues lie exactly on the imaginary axis. Apparently, this seems to
forbid the exponential stability obtained from Eq. (12) for the Lorentzian
distribution. The same mathematical conundrum was encountered and solved in
the Kuramoto model Strogatz and Mirollo (1991); Strogatz  _et al._ (1992);
Chiba and Nishikawa (2011), what suggests similar mechanisms are operating
here., we use the ansatz $\delta w(\eta,t)=b(\eta)e^{\lambda t}$ in (C2) and
obtain

$[\lambda+2w_{0}(\eta)i]b(\eta)e^{\lambda
t}=\frac{iJ}{2\pi}\int_{-\infty}^{\infty}[b(\eta)e^{\lambda
t}+b^{*}(\eta)e^{\lambda^{*}t}]g(\eta)d\eta.$

Solving the equation for $be^{\lambda t}$, we find

$b(\eta)e^{\lambda
t}=\frac{1}{2w_{0}(\eta)-i\lambda}\frac{J}{2\pi}\int_{-\infty}^{\infty}[b(\eta)e^{\lambda
t}+b^{*}(\eta)e^{\lambda^{*}t}]g(\eta)d\eta.$

After summing the complex conjugate at each side of the equation, we multiply
by $g(\eta)$ and integrate over $\eta$. This allows one to cancel out the
integrals $\int(be^{\lambda t}+b^{*}e^{\lambda^{*}t})gd\eta$. Then, imposing
self-consistency, we find

$1=\frac{J}{2\pi}\int_{-\infty}^{\infty}\left[\frac{1}{2w_{0}(\eta)-i\lambda}+\mbox{c.c.}\right]g(\eta)d\eta.$
(C3)

Local bifurcations are associated with marginal stability of the fixed points:
$\mathrm{Re}(\lambda)\to 0$. In correspondence with the results obtained for a
Lorentzian distribution in the main text, we expect bifurcations of saddle-
node type associated with $\lambda=0$. Therefore, from Eq. (C3) we find that
the loci of the saddle-node bifurcations are located at

$J_{SN}=\frac{2\pi}{\int_{-\infty}^{\infty}\frac{x_{0}(\eta)}{|w_{0}(\eta)|^{2}}g(\eta)d\eta}.$
(C4)

Since $x_{0}=\sqrt{\eta+Jr_{0}}$ for $\eta>-Jr_{0}$, and zero otherwise, we
can restrict the integration over $\eta$ to the range $(-Jr_{0},\infty)$. Now
multiplying Eq. (C4) by Eq. (C1) we get, after defining $\xi=J_{SN}r_{0}$, two
equations that permit to find the locus of the saddle-node bifurcations
systematically for arbitrary distributions of currents:

$J_{SN}=\frac{2\pi}{\int_{0}^{\infty}\eta^{-1/2}g(\eta-\xi)d\eta},$ (C5)

with $\xi$ obtained solving

$\xi=\frac{2\int_{0}^{\infty}\eta^{1/2}g(\eta-\xi)d\eta}{\int_{0}^{\infty}\eta^{-1/2}g(\eta-\xi)d\eta}.$
(C6)

Figure C1: (color online). Phase diagrams for (a) uniform and (b) Gaussian
distributions of currents $g(\eta)$. The cyan-shaded region represents the
bistable region, with the solid lines corresponding to saddle-node
bifurcations analytically obtained from Eqs. (C5) and (C6). Square symbols are
estimations of the bifurcations loci obtained from by direct numerical
simulations of $N=10^{4}$ QIF neurons. Triangle symbols indicate parameter
values used in numerical simulations of Figs. C2 and C3. Figure C2: Numerical
simulations of the excitatory network of QIF neurons, Eqs. (1), and (2). with
(a-d) uniform and (e-h) Gaussian distributions of currents. As in Fig. 2, at
time $t=0$, an external sinusoidal current $I(t)=I_{0}\sin(\omega t)$ —shown
in panels (d,h)— is applied to all neurons. Panels (a,e) show the time series
of the firing rate, and (b,f) the mean membrane potential of all cells. In
panels (c,g) we depicted the raster plots of 300 randomly chosen neurons.
Parameters correspond to the black triangle symbols in Fig. C1: (a-d)
$\bar{\eta}=-1$, $J=8$, $\gamma=1$, $I_{0}=1$; (e-h) $\bar{\eta}=-2$, $J=10$,
$\sigma=1$, $I_{0}=1.5$.

### Bistability region for uniform and Gaussian distributions

Equations (C5) and (C6) are particularly easy to solve for uniform
distributions of currents:

$g(\eta)=\begin{cases}\frac{1}{2\gamma}&\mbox{for}\,|\eta-\bar{\eta}|<\gamma\\\
0&\mbox{otherwise}\end{cases}$

After defining the rescaled parameters $\tilde{\eta}=\bar{\eta}/\gamma$ and
$\tilde{J}=J/\sqrt{\gamma}$, we find two branches of saddle-node bifurcations
emanating from the cusp point at $\tilde{\eta}=-1/3$:

$\displaystyle\tilde{J}_{SN}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{\sqrt{3\tilde{\eta}+3}}$
$\displaystyle\tilde{J}_{SN}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{\sqrt{\tilde{\eta}+1+2\sqrt{\frac{1}{3}+\tilde{\eta}^{2}}}-\sqrt{\tilde{\eta}-1+2\sqrt{\frac{1}{3}+\tilde{\eta}^{2}}}}$

These functions are plotted in Fig. S1(a). Numerical simulations using the
network of QIF neurons confirm the correctness of these boundaries (see the
square symbols in Fig. C1).

For Gaussian distributions,

$g(\eta)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(\eta-\bar{\eta})^{2}/(2\sigma^{2})},$

the solutions of (C5) and (C6) can be numerically found. The results are shown
in Fig. C1(b). Again, there is a perfect agreement between Eqs. (C5) and (C6)
—derived assuming the Lorentzian Ansatz— and the numerical estimations
obtained simulating the network of QIF neurons.

### Excitatory networks with external periodic currents

Firing rate equations (12) predict the existence of a stable focus in the
shaded region of Fig. 1. Trajectories attracted to this fixed point display
oscillations in the firing rate and mean membrane potential due to the
transient synchronous firing of the QIF neurons.

When an external periodic current is injected to all neurons in the network,
the spiral dynamics around the fixed point is responsible for the bursting
behavior observed in Fig. 2, as well as for the emergence of the macroscopic
chaos shown in Fig. 3. Here we investigate whether similar phenomena occur
when an external periodic current of the same frequency is injected in an
excitatory network with either uniform or Gaussian-distributed currents.

We introduced a sinusoidal forcing $I(t)=I_{0}\sin(\omega t)$ observing a
behavior qualitatively identical to the one reported in the main text with
Lorentzian $g(\eta)$. Under a low-frequency forcing the system displays
periodic bursting, see Fig. C2, provided the parameters are set inside the
bistable region (see the black triangles in Fig. C1). Furthermore, note that
the range of firing rates and mean membrane potentials in Figure C2 are
similar to those of Fig. 2 —for Lorentzian distributions of currents.

As for the simulations shown in Figure (3), we next increase the frequency of
the injected current up to $\omega=\pi$ to investigate the emergence of
macroscopic chaos in a network with Gaussian-distributed currents. Fig. C3
shows the emergence of an apparently chaotic state. The observed dynamics is
similar to that of Fig. 3, which suggests it is chaotic. This type of chaos
persists in the thermodynamic limit $N\to\infty$. On top of this, highly-
dimensional but weakly chaotic dynamics is probably also present due to
‘residual’ finite-size effects.

Figure C3: Chaotic state in a network of $N=10^{4}$ QIF neurons with Gaussian-
distributed currents. Network parameters correspond to the yellow-triangle
symbol in the phase diagram of Fig. S1(b). Like in Fig. 3, neurons receive a
common periodic current $I(t)=I_{0}\sin(\omega t)$ of frequency $\omega=\pi$.
Parameters: $\bar{\eta}=-1$, $J=8$, $\sigma=1$, $I_{0}=1.5$.

## Appendix D: Identical QIF neurons driven by independent noise terms

We compare now the results obtained above and in the main text for quenched
heterogeneity $g(\eta)$ with the results for identical neurons
$g(\eta)=\delta(\eta-\bar{\eta})$ driven by noise. Specifically, the inputs
currents are now taken as

$I_{j}=\bar{\eta}+Jr(t)+\xi_{j}(t),$ (D1)

where $\xi_{j}$ are independent white noise terms with expected values
$\left<\xi_{j}(t)\right>=0$, and
$\left<\xi_{j}(t)\xi_{k}(t^{\prime})\right>=2D\,\delta_{jk}\,\delta(t-t^{\prime})$.

Figure D1: (color online). Numerically obtained phase diagram (a) and time
series (b-d) for a network of identical QIF neurons driven by independent
white noise terms. The square symbols in panel (a) represent the boundaries
that enclose the region of bistability between a high and a low activity
state. The boundaries have been obtained using a network of $10^{4}$ neurons,
and by numerical continuation of a low activity solution (filled squares) and
a high activity solution (empty squares). Specifically, using a noise
intensity $D=1$ and a particular value of $J$, the system was initialized
either at $\bar{\eta}=-5$ or at $\bar{\eta}=0$ and, after $t=10$ time units,
parameter $\bar{\eta}$ was decreased/increased an amount $0.025/0.1$,
respectively. This continuation was made until the relative change of two
successive values of the averaged firing rate (time-averaged over the last
time unit for each parameter value) was larger than $50\%$. Other parameters
for the numerical simulations are the same as in all other figures, and are
described in the Material and Methods section. The triangle symbol indicates
the parameter value, $(\bar{\eta}/D^{2/3},J/D^{1/3})=(-5,15)$, corresponding
to the numerical simulations shown in black in the right panels —for the
simulations we used $D=1$. To facilitate the comparison with the FREs (12),
the orange curves show the time series of the FREs, using the same parameters
$(\bar{\eta},J)=(-5,15)$, but with $\Delta=1$

In the thermodynamic limit, the density $\rho(V,t)$ obeys the Fokker-Planck
equation:

$\partial_{t}\rho+\partial_{V}\left[(V^{2}+\bar{\eta}+Jr)\rho\right]=D\,\partial^{2}_{V}\rho.$
(D2)

The Lorentzian ansatz is not a solution of this equation, but we demonstrate
here that the phenomena observed with quenched heterogeneity arise also with
independent noise sources. This qualitative similarity at the macroscopic
level between quenched Lorentzian heterogeneity and Gaussian noise has been
noted in previous work Strogatz and Mirollo (1991).

A subtle point in Eq. (D2) is that its nondimensionalization entails a
different (compared to $\Delta$) rescaling with $D$: $\tilde{V}=V/D^{1/3}$,
$\tilde{\eta}=\eta/D^{2/3}$, $\tilde{J}=J/D^{1/3}$, $\tilde{t}=tD^{1/3}$
(implying $\tilde{r}=r/D^{1/3}$).

Numerical simulations reveal the existence of a region of bistability, see
Fig. D1(a), analogous to the one observed for quenched heterogeneity, cf.
Figs. 1(a) and C1. Obviously, true bistability only holds in the thermodynamic
limit, while what we observe are rather exceedingly long residence times close
to each fixed point (see figure caption for details).

Additionally, in order to investigate and compare the dynamical behavior of
the identical noise-driven neurons with that of the FREs (12), an external
current of intensity $I_{0}=3$ is injected to all neurons at time $t=0$, like
in Fig. 2. In Fig. D1(b,c) the time series of the firing rate and the mean
membrane potential clearly display damped oscillations after the injection of
the current, confirming the existence of a stable focus, exactly as observed
in the FREs (12). The existence of a stable focus reflects the presence of
transient spike synchrony in the network, as seen by the raster plot in Fig.
D1(d). Remarkably, the raster plot and the time series closely resemble those
of Fig. 2. The resemblance is not just qualitative, but rather there is a near
quantitative fit between the network of QIF neurons driven by Gaussian noise
and the FREs (12), which were derived assuming quenched Lorentzian
heterogeneity (orange curves).

## Appendix E: Model generalizations

To illustrate the potential of the LA for investigating more sophisticated
networks, here we provide the low-dimensional FREs corresponding to an
excitatory network of QIF neurons with independently distributed currents and
synaptic weights, and to two interacting populations of excitatory and
inhibitory QIF neurons with distributed currents.

### Excitatory population with heterogeneous currents and synaptic weights

As a first example, let us assume that both the currents $\eta$ and the
synaptic weights $J$ are distributed —this type of heterogeneity was also
considered in Buice and Chow (2013). The input currents read then 666Note that
considering stronger forms of heterogeneity for the synaptic weights may
generally not lead to Eq. (19).

$I_{j}=\eta_{j}+J_{j}r(t)+I(t).$

For simplicity, let us additionally assume that $\eta$ and $J$ are distributed
independently, with a joint distribution $p(\eta,J)=g(\eta)h(J)$. In the
simplest situation of Lorentzian $g(\eta)$ and $h(J)$

$g(\eta)=\frac{\Delta/\pi}{(\eta-\bar{\eta})^{2}+\Delta^{2}}\,;\qquad
h(J)=\frac{\Gamma/\pi}{(J-\bar{J})^{2}+\Gamma^{2}}.$

the problem is extremely simplified using the Lorentzian ansatz, which now
trivially reads:

$\rho(V|\eta,J,t)=\frac{1}{\pi}\frac{x(\eta,J,t)}{[V-y(\eta,J,t)]^{2}+x(\eta,J,t)^{2}}.$

The firing rate and mean membrane potential are determined only by the value
of $w\equiv x+iy$ at the poles in the lower half planes:

$\displaystyle r(t)+iv(t)$ $\displaystyle=$ $\displaystyle\iint
w(\eta,J,t)g(\eta)h(J)d\eta dJ$ $\displaystyle=$ $\displaystyle
w(\bar{\eta}-i\Delta,\bar{J}-i\Gamma,t).$

Finally, evaluating Eq. (10) in the main text at the poles
($\eta=\bar{\eta}-i\Delta$, $J=\bar{J}-i\Gamma$) we get the exact FREs:

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle\Delta/\pi+\Gamma
r/\pi+2rv,$ (E1a) $\displaystyle\dot{v}$ $\displaystyle=$ $\displaystyle
v^{2}+\bar{\eta}+\bar{J}r+I(t)-\pi^{2}r^{2}.$ (E1b)

These equations contain simply an extra term ‘$+\Gamma r/\pi$’ compared to
equations (12). Figure E1 shows the bistability boundaries obtained from (E1)
for different ratios of $\Gamma$ and $\Delta$, and $I(t)=0$. Note that the
region of bistability shifts to lower values of $\bar{\eta}/\Delta$ and to
higher values of $\bar{J}/\sqrt{\Delta}$ as the level of heterogeneity in the
synaptic coupling $\Gamma$ is increased.

Figure E1: (color online). Phase diagram for an excitatory population with
heterogeneous currents and synaptic weights obtained using the FREs (E1).

### Firing-rate equations for a pair of Excitatory-Inhibitory populations

The microscopic state of each population of QIF neurons is characterized by
the membrane potentials $\\{V_{j}^{(e,i)}\\}_{j=1,\ldots,N}$, which obey the
following ordinary differential equations:

$\dot{V}_{j}^{(e,i)}=(V_{j}^{(e,i)})^{2}+I_{j}^{(e,i)},\text{if
}V_{j}^{(e,i)}\geq V_{p},\text{then }V_{j}^{(e,i)}\leftarrow V_{r}.$

Here, $V_{j}^{(e,i)}$ represents the membrane potential of neuron $j$ in
either the excitatory $(e)$, or in the inhibitory population $(i)$. The
external currents for the excitatory and inhibitory populations are,
respectively:

$\displaystyle
I_{j}^{(e)}=\eta_{j}^{(e)}+J_{ee}s^{(e)}(t)-J_{ie}s^{(i)}(t)+I^{(e)}(t),$
$\displaystyle
I_{j}^{(i)}=\eta_{j}^{(i)}+J_{ei}s^{(e)}(t)-J_{ii}s^{(i)}(t)+I^{(i)}(t),$

where the synaptic weights are $J_{ee},J_{ii},J_{ie},J_{ei}$. Finally, the
mean synaptic activation for each population is

$s^{(e,i)}(t)=\frac{1}{N}\sum_{j=1}^{N}\sum_{k\backslash(t_{j}^{k})^{(e,i)}<t}\int_{-\infty}^{t}dt^{\prime}a_{\tau}(t-t^{\prime})\delta(t^{\prime}-(t_{j}^{k})^{(e,i)}).$

Here, $(t_{j}^{k})^{(e,i)}$ is the time of the $k$th spike of $j$th neuron in
either the excitatory $(e)$, or in the inhibitory population $(i)$.
Additionally, $\delta(t)$ is the Dirac delta function, and $a_{\tau}(t)$ is
the normalized synaptic activation caused by a single pre-synaptic spike with
time scale $\tau$, e.g. $a_{\tau}(t)=e^{-t/\tau}/\tau$.

It is straightforward to apply the LA and the method described in the main
text, to obtain the firing rate equations corresponding to the two-population
model. Considering the limit of infinitely fast synapses, $\tau\to 0$, we get
$s^{(e,i)}(t)=r^{(e,i)}(t)$. Finally, assuming the Lorentzian distributions of
currents for both populations:

$g_{e,i}(\eta)=\frac{1}{\pi}\frac{\Delta_{e,i}}{(\eta-\bar{\eta}_{e,i})^{2}+\Delta_{e,i}^{2}},$

we obtain the firing-rate equations:

$\displaystyle\dot{r}^{(e)}$ $\displaystyle=$
$\displaystyle\Delta_{e}/\pi+2r^{(e)}v^{(e)},$ $\displaystyle\dot{v}^{(e)}$
$\displaystyle=$
$\displaystyle(v^{(e)})^{2}+\bar{\eta}_{e}+J_{ee}r^{(e)}-J_{ie}r^{(i)}+I^{(e)}(t)-(\pi
r^{(e)})^{2},$ $\displaystyle\dot{r}^{(i)}$ $\displaystyle=$
$\displaystyle\Delta_{i}/\pi+2r^{(i)}v^{(i)},$ $\displaystyle\dot{v}^{(i)}$
$\displaystyle=$
$\displaystyle(v^{(i)})^{2}+\bar{\eta}_{i}+J_{ei}r^{(e)}-J_{ii}r^{(i)}+I^{(i)}(t)-(\pi
r^{(i)})^{2}.$