# Exact finite-dimensional description for networks of globally coupled
spiking neurons

Bastian Pietras<EMAIL_ADDRESS>Department of Information and
Communication Technologies, Universitat Pompeu Fabra, Tànger 122-140, 08018,
Barcelona, Spain  Rok Cestnik<EMAIL_ADDRESS>Department of Physics
and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25, 14476,
Potsdam-Golm, Germany  Arkady Pikovsky<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25,
14476, Potsdam-Golm, Germany

###### Abstract

We consider large networks of globally coupled spiking neurons and derive an
exact low-dimensional description of their collective dynamics in the
thermodynamic limit. Individual neurons are described by the Ermentrout-Kopell
canonical model that can be excitable or tonically spiking, and interact with
other neurons via pulses. Utilizing the equivalence of the quadratic
integrate-and-fire and the theta neuron formulations, we first derive the
dynamical equations in terms of the Kuramoto-Daido order parameters (Fourier
modes of the phase distribution) and relate them to two biophysically relevant
macroscopic observables, the firing rate and the mean voltage. For neurons
driven by Cauchy white noise or for Cauchy-Lorentz distributed input currents,
we adapt the results by Cestnik and Pikovsky [arXiv:2207.02302 (2022)] and
show that for arbitrary initial conditions the collective dynamics reduces to
six dimensions. We also prove that in this case the dynamics asymptotically
converges to a two-dimensional invariant manifold first discovered by Ott and
Antonsen. For identical, noise-free neurons, the dynamics reduces to three
dimensions, becoming equivalent to the Watanabe-Strogatz description. We
illustrate the exact six-dimensional dynamics outside the invariant manifold
by calculating nontrivial basins of different asymptotic regimes in a bistable
situation.

##  I Introduction

Synchronization phenomena in ensembles of coupled oscillators is an active
field of interdisciplinary research with numerous applications in physics,
engineering, and life sciences (see, e.g., books and reviews [1, 2, 3, 4]).
One important area of application is neuroscience, where synchronization of
neurons is crucial for understanding brain functioning [5, 6]. While for many
systems in physics and engineering the basic equations for the oscillator
dynamics can be formulated in the form of the Kuramoto model and its
modifications [3], neural models usually follow a different formulation that
takes into account specific properties of spiking neurons and their
interaction. Nevertheless, studies of collective effects in large neural
populations have profited enormously from the analogy of particular models of
neural dynamics to the Kuramoto model. Notably, a breakthrough in the
description of the Kuramoto-type dynamics by Ott and Antonsen [7] has recently
been transferred to the realm of spiking neuron networks [8, 9, 10, 11, 12],
see also the reviews [13, 14].

The main finding of Ott and Antonsen is the existence of an invariant two-
dimensional manifold, the so-called Ott-Antonsen (OA) manifold, corresponding
to a wrapped Cauchy distribution of the oscillators’ phases, which allows for
a formulation of exact closed equations of motion for the global order
parameter. In terms of the neural dynamics, these equations correspond to
closed equations for the parameters characterizing populations of neurons,
such as the firing rate and mean voltage. While the OA manifold is generally
assumed to be attractive for systems with quenched and/or noisy inhomogeneity,
and thus describes asymptotic in time attractors, generally initial conditions
lie outside the OA manifold and the corresponding transients are not captured
by the OA equations. Recently, two of us developed an exact six-dimensional
description for the evolution of Kuramoto-type oscillator populations outside
of the OA manifold [15]. The goal of this paper is to extend this approach to
networks of spiking neurons. We will demonstrate that with these exact six-
dimensional equations one can properly describe the transient dynamics outside
the OA manifold, and in particular find exact basins of attraction of
different asymptotic dynamical regimes (which lie on the OA manifold) for
arbitrary perturbations.

The paper is organized as follows. In Section II we introduce two basic
equivalent models of globally coupled spiking neurons, the theta-neurons and
the quadratic integrate-and-fire neurons, and derive equations for the order
parameters in the thermodynamic limit in presence of inhomogeneity and noise.
In Section III we discuss the mean field observables and the coupling terms.
Section IV contains our main findings. Here we derive the six-dimensional
reduction of the dynamics, and show how the initial states of the neurons
translate into initial conditions for the low-dimensional dynamics. We also
discuss the case of identical units (in the absence of noise and
heterogeneity). In Section V we discuss the relation to the OA theory, and in
particular we demonstrate how the stability of the OA manifold manifests
itself within our formalism. In Section VI we present several examples of
application of our approach to the dynamics off the OA manifold, including
numerical simulations of finite ensembles of spiking neurons. We discuss our
results in Section VII.

##  II Globally coupled spiking neurons

###  II.1 Theta-neuron and QIF formulations

In this paper we consider globally coupled spiking neurons described by the
Ermentrout-Kopell canonical model [16]. There are two equivalent formulations
of this model: one (theta-neuron, TN) uses a continuous phase-type variable
$\theta(t)$; the other (quadratic integrate-and-fire neuron, QIF) uses a
variable $v$ which is roughly interpreted as a membrane potential (for
mathematical simplicity one allows $v$ to attain infinite values, thus it is
not a real membrane potential of a neuron). We present both versions in
parallel, cf. [17, 8, 9, 10, 11, 12].

A QIF neuron is described by the equation

$\dot{v}=v^{2}-gv+I(t)\;,$ (1)

where $I(t)$ is the input current that together with the gap junction coupling
parameter $g$ determines the dynamics of the neuron. Because of the quadratic
nonlinearity, the voltage variable $v$ may reach infinity in a finite time, at
which point it is reset to $-\infty$, and this event is interpreted as a
spike.

To obtain the equation for a theta-neuron, one makes a transformation
$v(t)=\tan(\theta(t)/2)$, so that Eq. (1) takes on the form

$\dot{\theta}=1-\cos\theta-g\sin\theta+(1+\cos\theta)I(t)=\omega+\frac{he^{-i\theta}-h^{*}e^{i\theta}}{i}\;.$
(2)

Here, the spike occurs when the phase variable $\theta$ passes through $\pi$
with positive speed. Hence, the spike form can be modelled with a function of
the variable $\theta$ having a peak at $\theta\approx\pi$. We have introduced
parameters

$\omega(t)=1+I(t),\qquad h(t)=\frac{g+i[I(t)-1]}{2}\;,$ (3)

because they allow to interpret Eq. (2) as a Kuramoto-Sakaguchi model for
phase oscillators [18, 19, 3], with natural frequency $\omega$ and driving
field $h$.

The QIF dynamics (1) represents a discontinuous dynamical system that may
encompass real conceptual and mathematical difficulties due to the
instantaneous reset whenever a neurons crosses the threshold at infinity [20,
21]. For this reason, it may be advantageous to consider the TN dynamics (2),
which is bounded and smooth, circumvents the fire-and-reset discontinuity, and
can thus be treated within the theory of smooth dynamical systems.

###  II.2 Population of identical neurons in the thermodynamic limit

We now consider a population of identical neurons, which means that for all of
them the values of $g,I(t)$ are the same. In the thermodynamic limit of an
infinite number of units, a proper description is via the probability density.
In the TN formulation, the equation of probability conservation for this
density $P(\theta,t)$ reads

$\frac{\partial P}{\partial
t}+\frac{\partial}{\partial\theta}[P(\omega+\frac{1}{i}(he^{-i\theta}-h^{*}e^{i\theta}))]=0\;.$
(4)

In the QIF formulation, one has density $W(v,t)$ which is related to
$P(\theta,t)$ as

$W(v)=\frac{2P(2\arctan(v))}{1+v^{2}}\;,\quad
P(\theta)=\frac{W(\tan(\theta/2))}{1+\cos(\theta)}\;,$ (5)

(we have omitted the dependence on time for readability). This density obeys
the probability conservation equation

$\frac{\partial W}{\partial t}+\frac{\partial}{\partial
v}[W(v^{2}-gv+I)]=0\;,$

which should be equipped with the periodicity condition (outgoing flux at
“threshold” $v=\infty$ should be equal to the incoming flux at “reset”
$v=-\infty$).

The relation to the Kuramoto model mentioned above suggests the use of the
Kuramoto-Daido order parameters

$\displaystyle z_{n}(t)=\int_{0}^{2\pi}d\theta
P(\theta,t)e^{in\theta}=\int_{-\infty}^{\infty}dv\;W(v,t)\left(\frac{1+iv}{1-iv}\right)^{n},$
(6a) $\displaystyle
P(\theta,t)=\sum\limits_{n=-\infty}^{\infty}z_{n}(t)e^{-in\theta}\;$ (6b)

to describe the evolution, which now reduces to an infinite system

$\dot{z}_{n}=n[i\omega z_{n}+hz_{n-1}-h^{*}z_{n+1}]\;.$ (7)

This system is a Fourier-representation of Eq. (4). Since the probability
density is normalized, we have $z_{0}=1$.

The order parameter representation of the probability density deserves a more
detailed discussion. This representation is mostly straightforward in the TN
formulation, because there it is just the Fourier series representation of a
periodic function. Correspondingly, the order parameter description is
adequate if the probability density $P(\theta,t)$ is smooth enough [22]. For
the QIF model, the distribution density $W(v,t)$ is defined on the line
$-\infty<v<\infty$, and one may expect that a discrete series representation
would be insufficient. However, because we define the density $W$ via the
transformation (5), smoothness of $P$ implies a restriction on the possible
behavior of $W$ so that the series representation is valid for the QIF model
as well. This consideration shows that, in fact, the TN model is more adequate
for the description of spiking neurons. Indeed, the original derivation of the
TN model is based on the generic full equations of a Hodgin-Huxley-like
electrically excitable system [16]. By contrast, the QIF model has less
biophysical background and its interpretation should be based on the
transformation to the TN model.

###  II.3 Effect of Cauchy noise

Here we generalize the deterministic model above and consider spiking QIF
neurons (1) subject to additive independent, identically distributed Cauchy
noise $\xi(t)$ (cf. [23]); a particular case of this generalization is also
considered in [24]. The TN dynamics (2) now reads

$\dot{\theta}_{k}=1-\cos\theta_{k}-g\sin\theta_{k}+(1+\cos\theta_{k})\big{[}I(t)+\gamma\xi_{k}(t)\big{]}\;,$
(8)

where $\gamma$ is the noise intensity, and we assign index $k$ to neurons to
stress the independence of the noise terms $\xi_{k}$ and $\xi_{m}$ for $k\neq
m$. Then instead of Eq. (4) we obtain a generalized Fokker-Planck equation (in
the Stratonovich interpretation)

$\displaystyle\frac{\partial P}{\partial
t}+\frac{\partial}{\partial\theta}[P(1-\cos\theta-g$
$\displaystyle\sin\theta+(1+\cos\theta)I)]=$ (9) $\displaystyle=$
$\displaystyle\gamma\frac{\partial}{\partial|\theta|}(1+\cos\theta)P(\theta,t)\;,$

where in the Fourier space the operator $\frac{\partial}{\partial|\theta|}$
acts as $\frac{\partial}{\partial|\theta|}e^{ik\theta}=-|k|e^{ik\theta}$. If
we insert the Fourier series (6), then

$\displaystyle\frac{\partial}{\partial|\theta|}(1+\cos\theta)P(\theta,t)=$
$\displaystyle\sum\limits_{n=-\infty}^{\infty}z_{n}\Big{[}-|n|e^{-in\theta}-\frac{|n+1|}{2}e^{-i(n+1)\theta}-\frac{|n-1|}{2}e^{-i(n-1)\theta}\Big{]}\;.$

This leads to the following modification of Eqs. (7) for the order parameters
$z_{n},n\geq 1$,

$\begin{gathered}\dot{z}_{n}=n[i\omega
z_{n}+hz_{n-1}-h^{*}z_{n+1}-\gamma(z_{n}+\tfrac{1}{2}z_{n-1}+\tfrac{1}{2}z_{n+1})]\;.\end{gathered}$
(10)

###  II.4 Effect of Cauchy inhomogeneity

We further generalize the system under consideration by allowing for non-
identical neurons. We assume that the driving current $I(t)$ possesses an
additive quenched quantity $\eta_{k}$ so that instead of Eq. (8), we now
consider the TN dynamics

$\dot{\theta}_{k}=1-\cos\theta_{k}-g\sin\theta_{k}+(1+\cos\theta_{k})\big{[}I(t)+\gamma\xi_{k}(t)+\Delta\eta_{k}\big{]}\;,$
(11)

and the corresponding dynamics in the QIF formulation reads

$\dot{v}_{k}=v^{2}_{k}+I(t)+\gamma\xi_{k}(t)+\Delta\eta_{k}\;,$ (12)

with the same resetting as above. We now assume that the quantities $\eta_{k}$
(which are similar to natural frequencies in the Kuramoto model) are drawn
from a Cauchy-Lorentz distribution $f(\eta)$ with zero mean and unit half-
width,

$f(\eta)=\frac{1}{\pi(1+\eta^{2})}\;.$

The parameter $\Delta$ in Eqs. (11,12) determines the spread of “natural
frequencies”, i.e. the degree of heterogeneity among neurons.

Because the quantities $\eta_{k}$ are quenched (time-independent), we can
consider them as an additional parameter in the distribution of the variables
$\theta$, and write $P(\theta,t;\eta)$ instead of $P(\theta,t)$.
Correspondingly, one can introduce the order parameters $z_{n}(t;\eta)$ that
now depend also on the values of $\eta$. For globally coupled neurons, the
mean fields that are relevant then appear after averaging the order parameters
over the distribution of $\eta$:

$\displaystyle Z_{n}$
$\displaystyle=\int_{-\infty}^{\infty}d\eta\;f(\eta)z_{n}(t;\eta)$
$\displaystyle=\int_{-\infty}^{\infty}d\eta\,f(\eta)\int_{0}^{2\pi}d\theta\,P(\theta,t;\eta)e^{in\theta}\;.$

Next we follow the approach of Ott and Antonsen [7] and make an assumption
that the density $P(\theta,t;\eta)$ is an analytic function of the complex-
valued parameter $\eta$ in the upper half-plane and converges exponentially to
zero as $\text{Im}(\eta)\to\infty$. Then, the integral over $\eta$ can be
taken by virtue of the Cauchy residue theorem (using the pole $\eta=i$)

$\int_{-\infty}^{\infty}d\eta\,f(\eta)P(\theta,t;\eta)=P(\theta,t;i)\;,$

so that

$Z_{n}(t)=z_{n}(t;i)\;.$

Using Eqs. (10) together with $\omega=1+I(t)+\Delta\eta$ and
$h(t)=\frac{g+i(I(t)+\Delta\eta-1)}{2}$ evaluated at the pole $\eta=i$, we
finally obtain

$\begin{gathered}\dot{Z}_{n}=n\Big{[}i(1+I(t))Z_{n}-(\Delta+\gamma)Z_{n}+\hskip
34.14322pt\\\ +\frac{g+i(I(t)-1)-(\Delta+\gamma)}{2}Z_{n-1}-\hskip
14.22636pt\\\ -\frac{g-i(I(t)-1)+(\Delta+\gamma)}{2}Z_{n+1}\Big{]},\quad n\geq
1\;.\end{gathered}$ (13)

Noteworthy, both Cauchy noise and Cauchy inhomogeneity have the same effect on
the dynamics of the population, so that one can introduce the effective
inhomogeneity parameter $\Gamma=\Delta+\gamma$. We stress here, however, that
while the effect of noise in Eqs. (13) is unconditional, the effect of the
inhomogeneity is based on an additional assumption about analyticity of the
distribution density. This will be important for the interpretation of the
dynamical regimes below.

In another remark we would like to stress a special feature of the Cauchy-
distributed noise: because in Fourier space it acts proportionally to the mode
number, the resulting system (13) is homogeneous in mode number $n$, which is
essential for the theory below. This was first recognized in Refs. [23, 25].
In contradistinction, Gaussian noise acts proportionally to the square of the
mode number and thus destroys homogeneity of system (13).

##  III Mean field observables and coupling terms

In the description of populations of neurons, the two main observables are
typically the mean firing rate $R$ and mean voltage $V$. In networks of
globally coupled neurons, the synaptic input $I(t)$ depends on the recurrent
coupling and in some (but not all) cases becomes a function of $R$ and $V$.
Here, we relate the mean firing rate, mean voltage and synaptic input to the
order parameters.

###  III.1 Firing rate and mean voltage

Computational models of large networks of recurrently coupled spiking neurons
typically focus on a macroscopic observable that measures the mean rate at
which neurons emit spikes, the network firing rate

$R(t)=\frac{1}{N}\sum_{k=1}^{N}\sum_{j}\frac{1}{\tau_{r}}\int_{t-\tau_{r}}^{t}ds\;\delta(s-t_{k}^{(j)})\;,$
(14)

where the instant $t_{k}^{(j)}$ corresponds to the $j$-th spike of neuron $k$
(it happens when the corresponding variable $\theta_{k}$ crosses $\pi$), and
$\tau_{r}$ is a time window of spike events. Taking first the limit of
infinitely many neurons, $N\to\infty$, and then $\tau_{r}\to 0$, one obtains
the mean firing rate that, in terms of the probability density in Eq. (9), is
defined as the flux of probability density at $\theta=\pi$:

$R(t)=2P(\pi,t)=\frac{1}{\pi}\Big{[}1+\sum_{k=1}^{\infty}(-1)^{k}(Z_{k}+Z_{k}^{*})\Big{]}\;.$
(15)

The mean voltage $V$ is the population average of the membrane potential
variables $v$:

$V(t)=\langle
v\rangle=\left\langle\frac{\sin\theta}{1+\cos\theta}\right\rangle=i\sum_{k=1}^{\infty}(-1)^{k}(Z_{k}-Z_{k}^{*})\;.$
(16)

A combination of both $R$ and $V$ can be expressed as a simple alternating sum
of the moments:

$\pi R-iV=\sum\limits_{k=0}^{\infty}(-1)^{k}Z_{k}\;.$ (17)

###  III.2 Recurrent coupling

Next, we introduce two types of recurrent coupling in the population [26, 12]:
electrical coupling via gap junctions and chemical coupling via
excitatory/inhibitory synapses. Gap junction coupling of strength $g\geq 0$
amounts to including the term $gV$ to the driving current $I(t)$ of each
neuron [27]; for chemical interactions we include a term $\kappa S$, with
$\kappa>0$ $(\kappa<0)$ denoting excitatory (inhibitory) synaptic coupling.
The mean synaptic activity $S(t)$ satisfies the relaxation equation

$\tau_{s}\dot{S}=\langle\rho\rangle-S\;,$

where $\rho(\theta)$ is the spike pulse profile and $\tau_{s}$ a synaptic time
constant. In the limit of fast relaxation, i.e. for instantaneous interactions
$\tau_{s}\to 0$, one has
$S=\langle\rho\rangle=\int_{0}^{2\pi}d\theta\rho(\theta)P(\theta,t)$. Taken
together, the input current $I(t)$ consists of a common external input $I_{0}$
and the contributions from recurrent electrical and chemical coupling:

$I(t)=I_{0}+gV+\kappa S\;.$

In the following, we briefly present several possible choices of the pulse
profile $\rho(\theta)$. We assume that the pulse is localized around
$\theta=\pi$ (i.e. when the membrane potential diverges, $v\to\infty$), and
that the total area is normalized $\int_{0}^{2\pi}d\theta\rho(\theta)=2\pi$.
If we represent the pulse shape as a Fourier series

$\rho(\theta)=\sum_{k=-\infty}^{\infty}c_{k}e^{ik\theta},\quad c_{0}=1,$

then the average synaptic activity is represented via the order parameters as

$\langle\rho\rangle=1+\sum_{k=1}^{\infty}(c_{k}^{*}Z_{k}+c_{k}Z_{k}^{*})\;.$

#### Dirac delta-pulses.

Because spikes are rather narrow, in many situations a Dirac $\delta$-function
$\rho_{\delta}(\theta)=2\pi\delta(\theta-\pi)$ is adequate. The Fourier
coefficients are $c_{k}=(-1)^{|k|}$ and the synaptic activity reduces to the
mean firing rate (see Eq. (15))

$\langle\rho_{\delta}\rangle=1+\sum_{k=1}^{\infty}(-1)^{k}(Z_{k}+Z_{k}^{*})=\pi
R(t)\;.$ (18)

#### Ariaratnam-Strogatz (AS) pulse.

Ariaratnam and Strogatz [28] suggested a family of pulse profiles
$\rho_{AS}(\theta)=a_{m}(1-\cos\theta)^{m}$ with
$a_{m}=\frac{2^{m}(m!)^{2}}{(2m)!}$, where $m$ is an integer parameter. The AS
pulses are smooth, but in the limit $m\to\infty$ they converge to
$\delta$-pulses. The average synaptic activity is a finite sum of the order
parameters

$\langle\rho_{AS}\rangle=1+(m!)^{2}\sum_{k=1}^{m}(-1)^{k}\frac{Z_{k}+Z^{\ast}_{k}}{(m+k)!(m-k)!}\;.$
(19)

#### Rectified Poisson (RP) pulse.

Gallego et al. [29] suggested the following pulse form

$\rho_{RP}(\theta)=\frac{(1-r)(1-\cos\theta)}{1+2r\cos\theta+r^{2}}\;,$

where the “sharpness” parameter $r\in(-1,1)$ determines the width of the
pulse: For $r=-1$, $\rho_{RP}=1$ is flat; for $r=0$,
$\rho_{RP}(\theta)=1-\cos\theta$ coincides with the AS pulse with $m=1$; and
in the limit $r\to 1$, $\rho_{RP}(\theta)$ becomes a $\delta$-pulse. The
Fourier coefficients are $c_{k}=\frac{1+r}{2r}(-r)^{|k|}$ and the average
activity is the infinite series

$\langle\rho_{RP}\rangle=1+\frac{1+r}{2r}\sum_{k=1}^{\infty}(-r)^{k}\big{(}Z_{k}+Z^{\ast}_{k}\big{)}\;.$
(20)

The different pulse shapes above are symmetric about $\theta=\pi$ and satisfy
$\rho(0)=0$. We remark that the RP pulses of width $r$ can be further
generalized to account for non-symmetric pulses when considering Fourier
coefficients $c_{k}=c_{-k}^{*}$ of the form
$c_{k}=ae^{i(\varphi+k\psi)}r^{k},k>0,$ with additional parameters $a>0$ and
$\varphi,\psi\in[0,2\pi)$. In all these cases, it is possible to represent the
relevant observables and mean fields governing the dynamics of the neural
population via the order parameters $Z_{k}$. Such representations can directly
be used for numerical simulations of the ensemble of neurons in the
thermodynamic limit, with a proper truncation of the infinite series.

##  IV Finite-dimensional reduction

###  IV.1 Reduction in terms of Kuramoto order parameters

In this section we propose that the infinite system of order parameter
dynamics (13) can be reduced to three complex equations (and a constant
function, which is an integral of motion). Our derivation here directly
follows the corresponding derivation for the Kuramoto problem in [15], so we
omit some technical details.

First, we rewrite Eqs. (13) in a more compact form

$\dot{Z}_{n}=n\big{[}(i\omega-\Gamma)Z_{n}+(h-\tfrac{\Gamma}{2})Z_{n-1}-(h^{*}+\tfrac{\Gamma}{2})Z_{n+1}\big{]},\;n\geq
1.$ (21)

We introduce the complex-valued exponential generating function (EGF)
$\mathbf{Z}(k,t)=\sum_{n=0}^{\infty}Z_{n}(t)\frac{k^{n}}{n!}$, which obeys the
partial differential equation (a prime denotes derivative with respect to
$k$):

$\dot{\mathbf{Z}}=k(i\omega-\Gamma)\mathbf{Z}^{\prime}+k(h-\tfrac{\Gamma}{2})\mathbf{Z}-k(h^{*}+\tfrac{\Gamma}{2})\mathbf{Z}^{\prime\prime}\;.$

With the Ansatz $\mathbf{Z}(k,t)=e^{kQ(t)}\mathbf{B}(k,t)$ assuming that $Q$
obeys $\dot{Q}=i\omega Q+h-h^{*}Q^{2}-\frac{\Gamma}{2}(1+Q)^{2}$, we find that
the generating function
$\mathbf{B}(k,t)=\sum_{n=0}^{\infty}\beta_{n}(t)\frac{k^{n}}{n!}$ obeys the
dynamics

$\dot{\mathbf{B}}=k\left[i\omega-\Gamma-2(h^{*}+\tfrac{\Gamma}{2})Q\right]\mathbf{B}^{\prime}-k(h^{*}+\tfrac{\Gamma}{2})\mathbf{B}^{\prime\prime}$

and the equations for the new variables $\beta_{n}$ read

$\tfrac{1}{n}\dot{\beta}_{n}=\left[i\omega-\Gamma-2(h^{*}+\tfrac{\Gamma}{2})Q\right]\beta_{n}-(h^{*}+\tfrac{\Gamma}{2})\beta_{n+1},\;n\geq
1.$ (22)

Due to normalization of the exponential generating function
$\mathbf{Z}(0,t)=1$, we have $\beta_{0}=1$.

We now introduce two new complex dynamical variables $y$ and $s$. As we will
show below, the set of variables $\beta_{n}$ can fully be represented through
these variables $y,s$ and the constants of motion to be defined below. The
collective dynamics of the network of spiking neurons is exactly described by
the following dynamical equations for $Q,y,s$,

$\displaystyle\dot{Q}$ $\displaystyle=i\omega
Q+h-h^{\ast}Q^{2}-\tfrac{\Gamma}{2}(1+Q)^{2}\;,$ (23a) $\displaystyle\dot{y}$
$\displaystyle=\big{[}i\omega-2h^{\ast}Q\big{]}y-\Gamma(1+Q)y\;,$ (23b)
$\displaystyle\dot{s}$ $\displaystyle=h^{\ast}y+\tfrac{\Gamma}{2}y\;.$ (23c)

To reveal the connection between $\\{\beta_{n}\\}$ and $y,s$, we introduce an
additional set of variables $\\{\alpha_{n}\\}$ via
$\beta_{n}=y^{n}\alpha_{n}$. In combination with Eqs. (22) and (23b), we find
the dynamics of $\alpha_{n}$

$\tfrac{1}{n}\dot{\alpha}_{n}=-(h^{*}+\tfrac{\Gamma}{2})y\alpha_{n+1},\quad
n\geq 1\;.$ (24)

In terms of the corresponding ordinary generating function (OGF)
$\mathcal{A}(k,t)=\sum_{n=1}^{\infty}\alpha_{n}k^{n}$, the dynamics (24) can
be represented as

$\dot{\mathcal{A}}=-(h^{*}+\tfrac{\Gamma}{2})y\left[\mathcal{A}^{\prime}-k^{-1}\mathcal{A}\right]\;.$
(25)

Finally, we introduce another OGF that is related to $\mathcal{A}$ via the
dynamical variable $s(t)$:

$\mathcal{M}(k)=\frac{k}{k+s(t)}\mathcal{A}\big{(}k+s(t),t\big{)}\;.$

A direct calculation of the time derivative of $\mathcal{M}$ yields
$\dot{\mathcal{M}}(k)=0$, which means that this function is an integral of
motion. Equivalently, using the representation
$\mathcal{M}(k)=\sum_{n=1}^{\infty}\mu_{n}k^{n}$, we observe that the system
possesses an infinite set of integrals $\mu_{n}$. To conclude the derivation,
we present the relations between the different variables (see [15] for an
explicit derivation):

$\displaystyle\mu_{n}$
$\displaystyle=\sum\limits_{m=n}^{\infty}\binom{m-1}{n-1}\alpha_{m}(t)\big{[}s(t)\big{]}^{m-n}\;,$
(26a) $\displaystyle\alpha_{n}$
$\displaystyle=\sum\limits_{m=n}^{\infty}\binom{m-1}{n-1}\mu_{m}\big{[}-s(t)\big{]}^{m-n}\;,$
(26b) $\displaystyle Z_{n}$
$\displaystyle=\sum\limits_{m=0}^{n}\binom{n}{m}\beta_{m}(t)\big{[}Q(t)\big{]}^{n-m}$
(26c)
$\displaystyle=Q^{n}-\sum\limits_{m=1}^{n}\binom{n}{m}Q^{n-m}y^{m}\sum\limits_{d=0}^{m-1}\frac{s^{d-m}}{d!}\mathcal{M}^{(d)}(-s)\;,$

where $\mathcal{M}^{(d)}$ denotes the $d^{\text{th}}$ derivative of
$\mathcal{M}$ with respect to $k$. The Kuramoto order parameter, i.e. the
first circular moment, expresses as:

$Z_{1}=Q-y\frac{\mathcal{M}(-s)}{s}\;.$ (27)

###  IV.2 Reduction in terms of voltage and firing rate

As we discussed in Section III, macroscopic observables of interest for
populations of neurons are the mean firing rate $R$ and the mean voltage $V$.
It is instructive to represent them in terms of the newly introduced dynamical
variables $Q,y,s$. Inserting Eq. (26c) into the expansions (15) and (16) for
$R$ and $V$, respectively, and using the definitions of the OGF
$\mathcal{A,M}$, after straightforward algebra we obtain

$\displaystyle R$
$\displaystyle=\frac{1}{\pi}\text{Re}\left[\frac{1-Q}{1+Q}+\frac{2y}{y(1+Q)+s(1+Q)^{2}}\mathcal{M}\left(-\frac{y+s(1+Q)}{1+Q}\right)\right]\;,$
(28a) $\displaystyle V$
$\displaystyle=-\text{Im}\left[\frac{1-Q}{1+Q}+\frac{2y}{y(1+Q)+s(1+Q)^{2}}\mathcal{M}\left(-\frac{y+s(1+Q)}{1+Q}\right)\right]\;.$
(28b)

Expressions (28) suggest a transformation from the variables $Q,y,s$ to new
dynamical variables $\Phi,\lambda,\sigma$ according to

$\Phi=\frac{1-Q}{1+Q},\quad\lambda=\frac{2y}{(1+Q)^{2}},\quad\sigma=\frac{y}{1+Q}+s\;.$
(29)

The inverse transformation reads

$Q=\frac{1-\Phi}{1+\Phi},\quad y=\frac{2\lambda}{(1+\Phi)^{2}},\quad
s=\sigma-\frac{\lambda}{1+\Phi}\;.$ (30)

When substituting Eq. (3) for $\omega,h$, we obtain the dynamical equations
for the new variables

$\displaystyle\dot{\Phi}$ $\displaystyle=i\Phi^{2}-g\Phi-iI(t)+\Gamma\;,$
(31a) $\displaystyle\dot{\lambda}$ $\displaystyle=2i\Phi\lambda-g\lambda\;,$
(31b) $\displaystyle\dot{\sigma}$ $\displaystyle=i\lambda\;.$ (31c)

The main advantage is that the mean firing rate and the mean voltage in terms
of these new variables simplify as

$\pi R-iV=\Phi+\lambda\frac{\mathcal{M}(-\sigma)}{\sigma}\;.$ (32)

At this point some comments are due: First, note the resemblance between Eqs.
(32) and (27), which may be due to the transformal mapping between $\Phi$ and
$Q$. Second, there is also a more apparent resemblance of Eq. (31a) to the
initial QIF dynamics (1). Third, Eqs. (31) are asymmetrically coupled: the
dynamics (31a) for $\Phi$ does not depend on $\lambda$ and $\sigma$
explicitly, but may implicitly depend on them through the input current $I$,
e.g., if $I=I(R,V)$ depends on mean-field terms such as the firing rate $R$ or
the mean voltage $V$. In case that the current $I$ is not a function of the
state variables $\lambda$, $\sigma$, then Eqs. (31) is a skew system and the
first equation (31a) acts as a two-dimensional driving to the other two
equations.

###  IV.3 Recurrent coupling in terms of dynamical variables

For a self-consistent description of the collective dynamics, we also need to
express the recurrent coupling via gap junctions and/or via the previously
introduced synaptic pulses in terms of the newly introduced variables
$\Phi,\lambda,\sigma$. As the gap junction coupling term is proportional to
the mean voltage $V$, we can here use Eq. (32). Likewise, for synaptic
interactions via Dirac $\delta$-pulses $\rho_{\delta}(\theta)$ according to
Eq. (18), we can also use Eq. (32) because the average synaptic activity
$\langle\rho_{\delta}\rangle$ reduces to the mean firing rate $R$.

In general, however, the synaptic coupling depends on the assumed shape of the
pulse and one obtains complex expressions for $\langle\rho\rangle$, e.g., for
the RP pulse (20) we have

$\displaystyle\langle\rho_{RP}\rangle$
$\displaystyle=\text{Re}\left[\frac{1-Q}{1+rQ}+\frac{(1+r)y}{ry(1+rQ)+s(1+rQ)^{2}}\mathcal{M}\left(-\frac{ry+s(1+rQ)}{1+rQ}\right)\right]$
$\displaystyle=\text{Re}\left[\frac{\Phi}{u}+\frac{(1+r)\lambda}{2\sigma
u^{2}-(1-r)\lambda
u}\mathcal{M}\left(-\sigma+\frac{(1-r)\lambda}{u}\right)\right]\;,$

where we introduced the auxiliary variable
$u=\frac{1}{2}(1+r+(1-r)\Phi)=\frac{1+rQ}{1+Q}$ with $u\big{|}_{r=1}=1$.

For the AS pulse, $\langle\rho_{AS}\rangle$ is represented by a finite sum of
order parameters (19). Here one has to use Eq. (26c) expressing these order
parameters via dynamical variables $Q,y,s$. If needed, the transformation (30)
to variables $\Phi,\lambda,\sigma$ can be used.

###  IV.4 Initial conditions

Given an initial distribution $P(\theta,t=0)$ of phases $\theta(0)$, or
$W(v,t=0)$ of membrane potentials $v(0)$, we ought to initialize the dynamical
system $(Q,y,s)$ accordingly. Note however, that the set of dynamical
variables $Q,y,s$ is underdetermined, so that there is a certain freedom in
choosing the initial values $Q(0),y(0),s(0)$. We adopt the approach of [15]
and define the initial values as

$Q(0)=s(0)=0,\quad y(0)=1.$ (33)

In this case, as it follows from relations (26) in Section IV.1,

$\mu_{n}=\alpha_{n}(0)=\beta_{n}(0)=Z_{n}(0)\;.$

The constant OGF $\mathcal{M}(k)$ is thus directly related to the initial
values of the order parameters

$\mathcal{M}(k)=\sum_{n=1}^{\infty}Z_{n}(0)k^{n}\;.$

Substituting here the expression of the order parameters via the initial
distribution density of the phase variables $P(\theta,0)$, we obtain

$\begin{gathered}\mathcal{M}(k)=\int\limits_{0}^{2\pi}d\theta\;P(\theta,0)\frac{ke^{i\theta}}{1-ke^{i\theta}}=\\\
=\int\limits_{-\infty}^{\infty}dv\;W(v,0)\frac{k(1+iv)}{1-k-iv(1+k)}\;.\end{gathered}$
(34)

In terms of the observables $\Phi,\lambda,\sigma$ the initial conditions are

$\Phi(0)=\sigma(0)=1,\quad\lambda(0)=2\;.$ (35)

Below we present several cases where an expression for $\mathcal{M}(k)$ is
relatively simple.

  * •

A delta-distribution of voltages $W(v,0)=\delta(v-v_{0})$ corresponds to
$P(\theta,0)=\delta(\theta-\theta_{0})$ with $\theta_{0}=2\arctan(v_{0})$,
leading to $Z_{n}=e^{in\theta_{0}}$ and thus
$\mathcal{M}(k)=\frac{e^{i\theta_{0}}k}{1-e^{i\theta_{0}}k}=k\big{(}\frac{1+iv_{0}}{1-iv_{0}}-k\big{)}^{-1}.$

  * •

A Cauchy-Lorentz distribution of voltages

$W(v,0)=\frac{1}{\pi}\frac{x}{(v-v_{0})^{2}+x^{2}}\text{ for some
}x\in\mathbb{R}_{>0}$ (36)

corresponds to a wrapped Cauchy distribution of the phase variables
$P(\theta,0)=\frac{1}{2\pi}\frac{1-|\mu|^{2}}{|1-\mu e^{-i\theta}|^{2}}$ where
$\mu=\frac{1-\pi x+iv_{0}}{1+\pi x-iv_{0}}\in\mathbb{C}$. This leads to
$\mathcal{M}(k)=\frac{\mu k}{1-\mu k}$. The initial order parameters are
powers of the parameter, $Z_{n}(0)=\mu_{n}=\mu^{n}$.

A skewed generalization of the Cauchy distribution can also easily be
represented with the Kato-Jones [30] distribution of phases, which yields
$\mathcal{M}(k)=c\frac{\mu k}{1-\mu k}$ where $c\in\mathbb{C}$.

  * •

A uniform distribution of voltages in some interval $[v_{0}-d,v_{0}+d]$ with
$v_{0}\in\mathbb{R},d>0$,

$W(v,0)=\begin{cases}\frac{1}{2d}\quad&\text{if }v\in[v_{0}-d,v_{0}+d],\\\
0,&\text{otherwise},\end{cases}$ (37)

corresponds to $P(\theta,0)=1/(2d)/(1+\cos\theta)$ if
$-\pi<2\arctan(v_{0}-d)\leq\theta\leq 2\arctan(v_{0}+d)<\pi$ and $0$
otherwise, where we used that $\arctan(v)$ is bijective on $(-\pi,\pi)$. The
OGF is

$\displaystyle\mathcal{M}(k)=\frac{-1}{d}\frac{k}{(1+k)^{2}}\Big{[}(1+k)d+$
(38) $\displaystyle\hskip
11.38092pt+i\log\left(\frac{(1+k)(v_{0}-d+i)-2ik}{(1+k)(v_{0}+d+i)-2ik}\right)\Big{]}.$

We notice also that a weighted sum of different initial distributions allows
for a proper representation of the constant OGF $\mathcal{M}(k)$. Because of
linearity of representation (34), if the initial density of the phase variable
is a sum of “elementary” densities, then the OGF is a sum of “elementary”
OGFs:

$P(\theta,0)=\sum_{m}w_{m}P_{m}(\theta,0),\quad\mathcal{M}(k)=\sum_{m}w_{m}\mathcal{M}_{m}(k)$
(39)

for weights $w_{m}\geq 0$ with $\sum_{m}w_{m}=1$. We will explore this
property in the next Section IV.5 and in the Section VI below, where we
discuss the dynamics of a population of neurons following the resetting of a
fraction of neurons to a new state.

Finally, we mention that the function $\mathcal{M}(k)$ can be approximated for
an arbitrary distribution $W(v,0)$, by first expressing it in the TN
formulation $P(\theta,0)$ via (5), and then approximating it with a finite
Fourier series:
$P(\theta)\approx\frac{1}{2\pi}\big{[}1+\sum\limits_{n=1}^{N}a_{n}e^{-in\theta}+\text{c.c.}\big{]}$,
yielding a polynomial $\mathcal{M}(k)\approx\sum\limits_{n=1}^{N}a_{n}k^{n}$.

###  IV.5 Identical units without noise

In the case of identical units without noise, $\Gamma=0$, the following
simplification follows from Eq. (35) in [15]:

$\lambda=(\Phi+\Phi^{*})\sigma\;.$ (40)

The dynamics (31b) for $\lambda$ thus becomes redundant and Eq. (31c) for
$\sigma$ reduces to
$\dot{\sigma}=i(\Phi+\Phi^{*})\sigma=2i\text{Re}[\Phi]\sigma$. Given the
initial value $\sigma(0)=1$, we thus have $\sigma(t)=\exp[i\zeta(t)]$ and the
argument $\zeta\equiv\arg(\sigma)$ follows the real-valued dynamics
$\dot{\zeta}=2\text{Re}[\Phi]$. Hence, for $\Gamma=0$ the collective dynamics
is 3-dimensional. Note that the relation (32) connecting the firing rate and
mean voltage remains valid as well as the different choices of initial
conditions described in Section IV.4. The full three-dimensional dynamics can
be summarized as:

$\displaystyle\dot{\Phi}$ $\displaystyle=i\Phi^{2}-g\Phi-iI(t)\;,$ (41a)
$\displaystyle\dot{\zeta}$ $\displaystyle=2\text{Re}[\Phi]\;,$ (41b)

with initial conditions $\Phi(0)=1$, $\zeta(0)=0$; and the firing rate $R$ and
mean voltage $V$ can be expressed via Eq. (32) as:

$\pi R-iV=\Phi+(\Phi+\Phi^{*})\mathcal{M}(-e^{i\zeta})\;.$ (42)

This essentially boils down to the Watanabe-Strogatz (WS) description [31, 32]
of a population of identical units, here translated for QIF neurons. The
constant function $\mathcal{M}(k)$ contains the information of the WS
constants of motion $\psi_{j}=\theta_{j}(0)$, which are, for this simple
choice of initial conditions (35), just the initial values of the TN phases
and relate to the initial voltages via $v_{j}(0)=\tan(\psi_{j}/2)$. The
function $\mathcal{M}(k)$ is in general expressed as:

$\mathcal{M}(k)=\sum\limits_{n=1}^{\infty}\left\langle{e^{in\psi_{j}}}\right\rangle
k^{n}=\sum\limits_{n=1}^{\infty}\left\langle{\left(\frac{1+iv_{j}(0)}{1-iv_{j}(0)}\right)^{n}}\right\rangle
k^{n}\;.$

This description outlined above is also applicable for finite networks:
Consider an ensemble of $N$ identical neurons
$\\{\theta_{1},\theta_{2},...,\theta_{N}\\}$ (or $\\{v_{1},v_{2},...,v_{N}\\}$
in the QIF representation), whose distribution density can be described as a
sum of delta functions:
$P(\theta)=\frac{1}{N}\sum\limits_{n=1}^{N}\delta(\theta-\theta_{n})$ (or
$W(v)=\frac{1}{N}\sum\limits_{n=1}^{N}\delta(v-v_{n})$). For each unit
$\theta_{n}$ (or $v_{n}$), we can determine the constant function
$\mathcal{M}_{n}(k)$ corresponding to the respective delta function, that is,
$\mathcal{M}_{n}(k)=\frac{e^{i\theta_{n}(0)}k}{1-e^{i\theta_{n}(0)}k}=k\big{(}\frac{1-iv_{n}(0)}{1+iv_{n}(0)}-k\big{)}^{-1}$,
$n=1,\dots,N$. Using Eq. (39), we thus obtain the full function
$\mathcal{M}(k)=\frac{1}{N}\sum\limits_{n=1}^{N}\mathcal{M}_{n}(k)$. In this
way, ensembles with an arbitrary number of QIF neurons can be integrated
exactly in just 3 dimensions as their collective dynamics is restricted to the
three-dimensional manifold with coordinates $|\Phi|,\arg(\Phi),\zeta$. We
stress again that this three-dimensional description is only valid for
identical neurons without noise, $\Gamma=0$, i.e. without dissipation.

##  V Ott-Antonsen manifold and its stability

###  V.1 Ott-Antonsen and Lorentzian manifolds

From Eqs. (23) it is obvious that $y=0$ is a solution. The manifold
$\\{y=0\\}$ is invariant, and on it the variable $s$ becomes irrelevant
because, as it follows from Eq. (26c), the order parameters are just powers of
the variable $Q$: $Z_{n}=Q^{n}$. This manifold is thus completely described by
the complex variable $Q$, and hence is two-dimensional. It has first been
identified by Ott and Antonsen [7], who subsequently demonstrated that, in the
case of heterogeneous oscillators $\omega_{j}\neq\omega_{k}$ with common
forcing $h_{j}=h$, it is attractive in the time-asymptotic limit [33, 34]. Due
to the relation $Z_{n}=Q^{n}$, the distribution $P(v,t)$ of phases $\theta$ on
the so-called OA manifold follows a wrapped Cauchy distribution (a.k.a.
Poisson kernel).

In terms of the variables $\Phi,\lambda,\sigma$, the OA manifold corresponds
to $\\{\lambda=0\\}$ and on it the dynamics are fully captured by the variable
$\Phi=\pi R-iV$. Noteworthy, for $\lambda=0$ the distribution of voltages $v$
is a Cauchy-Lorentz one, thus in this context one calls the OA manifold also
the Lorentzian manifold [11]. The dynamics of neural populations on the OA, or
Lorentzian, manifold has been explored in a variety of setups, building
primarily on the original results in Refs. [8, 9, 10, 11, 12]. Time-asymptotic
attractiveness of the Lorentzian manifold was later proven in [35].

###  V.2 Global stability of the reduced OA manifold

Here we argue that the reduced OA manifold for TN, or the Lorentzian manifold
for QIF neurons, is asymptotically stable for $\Gamma>0$, i.e. in presence of
an inhomogeneity of input currents $\eta$ or of Cauchy white noise $\xi$. In
fact, we show that asymptotically $y\to 0$ so that the OA manifold is
attractive. Our starting point are Eqs. (23a),(23b). We introduce two real
variables $X=1-|Q|^{2}$ and $Y=|y|^{2}$, and express their dynamics as

$\displaystyle\frac{\dot{X}}{X}$
$\displaystyle=-(h^{*}Q+hQ^{*})+2\Gamma\frac{1-X}{X}+\Gamma(Q+Q^{*})\frac{2-X}{2X}\;,$
(43a) $\displaystyle\frac{\dot{Y}}{Y}$
$\displaystyle=-(h^{*}Q+hQ^{*})-\Gamma(1+\frac{Q+Q^{*}}{2})\;.$ (43b)

First, we argue that $X$ cannot become negative. Indeed, because at $X=0$,
$|Q|=1$, $2+Q+Q^{*}\geq 0$ and therefore $X$ cannot become negative. Moreover,
$X$ cannot vanish because for this one needs the state with $Q=-1$ to be
steady, but at this value of $Q$ from Eq. (23a) we have
$\dot{Q}\big{|}_{Q=-1}=-2i$ independently of the coupling terms $g,I$. This
corresponds to the fact that in the original formulation (2) at $\theta=\pi$
one has $\dot{\theta}=2$ independently of the forcing terms $g,I$.

Next, we combine these two equations in one:

$\frac{\dot{Y}}{Y}=\frac{\dot{X}}{X}-\Gamma\frac{|1+Q|^{2}}{1-|Q|^{2}}\;.$
(44)

We integrate this equation to obtain

$Y(t)=Y(0)\frac{X(t)}{X(0)}\exp[-\Gamma\int\limits_{0}^{t}\frac{|1+Q(t^{\prime})|^{2}}{1-|Q(t^{\prime})|^{2}}\,dt^{\prime}]\;.$
(45)

If $|Q(t)+1|$ is bounded from zero, then also
$\frac{|1+Q(t)|^{2}}{1-|Q(t)|^{2}}$ has a lower positive bound, and so the
integral in Eq. (45) tends to zero exponentially for $\Gamma>0$. On the other
hand, if $Q\to-1$, then $X$ vanishes and in this case $Y$ vanishes as well.
Thus in all the cases $Y(t)$ eventually vanishes, which means convergence of
arbitrary initial states to the OA manifold. Furthermore, $Y\to 0$ entails
$y\to 0$, which in turn entails $\lambda\to 0$ due to Eq. (29), thus proving
global stability also of the Lorentzian manifold described by $\Phi$.

The convergence $Y\to 0$ is ensured for $\Gamma>0$, but for vanishing noise
and/or inhomogeneity $\Gamma=0$ we obtain $Y(t)=X(t)\frac{Y(0)}{X(0)}$. Thus,
$Y$ can vanish only if eventually $X\to 0$, which corresponds to a delta-
function distribution with $|Q|=1$, i.e. to full synchrony of the neurons. In
all other situations the OA manifold is not attracting.

##  VI Dynamics off the OA manifold

As we have discussed in the previous Section V, asymptotically the six-
dimensional dynamics reduces to a two-dimensional dynamics on the OA manifold.
Thus, the exact six-dimensional evolution derived above is relevant for
transient processes only, as the attractors themselves lie on the OA manifold.
Importantly, the full six-dimensional dynamics is needed to faithfully
determine the basins of attraction of different asymptotic regimes. The two-
dimensional OA approach is restricted to specific initial conditions that
already lie on the invariant OA manifold (wrapped Cauchy distribution for TN
phases, Cauchy-Lorentz distribution for QIF voltages). By contrast, our
approach captures the exact six-dimensional evolution from general initial
phase/voltage distributions, which are incorporated in the corresponding
$\mathcal{M}$-function as described in Section IV.4. In the following, we
present three examples that underscore the power of our exact finite-
dimensional reduction.

Figure 1: Comparison between network dynamics of $N=20\,000$ neurons, the
exact six-dimensional collective dynamics, and the two-dimensional OA theory.
Collective dynamics start from a uniform distribution of initial voltages
$v_{j}(0)\in[v_{0}-d,v_{0}+d]$ with $d=1$ and (a) $v_{0}=0.75$ or (b)
$v_{0}=0.25$. Top: raster plot of $1000$ randomly chosen spiking neurons with
dynamics (8). Perfect agreement between network firing rate $R(t)$ (middle) /
mean voltage $V(t)$ (bottom) and exact low-dimensional theory Eqs. (31),
whereas the OA dynamics run into a different attractor. Network firing rate
and voltage are smoothed with a rectangular filter of width $0.05$. Other
parameters: noise strength $\gamma=1/4$, heterogeneity $\Delta=0$, i.e.
$\Gamma=1/4$; coupling strengths $\kappa=15\sqrt{\Gamma}/\pi$, $g=0$; input
current $I_{0}=-4\Gamma$.

###  VI.1 Complex initial transients of collective dynamics

As a first litmus test, we compare our finite-dimensional reduction Eqs. (31)
to numerical simulations of large networks of spiking neurons, see Fig. 1. For
the network simulations, we consider $N=20\,000$ neurons driven by Cauchy
white noise, whose dynamics is given by Eq. (8). The neurons are all-to-all
coupled via instantaneous Dirac $\delta$-pulses (Eq. (18)) of strength
$\kappa>0$, and $g=0$. The initial voltages $v_{j}(0)$ are uniformly
distributed in $[v_{0}-1,v_{0}+1]$, so that the initial firing rate $R(0)=0$
and mean voltage $V(0)=v_{0}$. The population firing rate $R(t)$ can be
computed according to Eq. (14), or in the TN framework as
$R(t)=\frac{1}{N\;dt}\sum_{j=1}^{N}1_{(\pi-2dt,\pi]}\big{(}\theta_{j}(t)\big{)}$,
where $1_{A}(\theta)$ is the indicator function and the interval
$A=(\pi-2dt,\pi]$ ensures that we count all neurons that cross the firing
threshold $\theta=\pi$ within the next integration step $dt$, and thus elicit
a spike before time $t+dt$. The mean voltage $V(t)$ is computed as the average
over the neurons’ voltages, or in the TN framework as
$V\approx\langle\frac{\sin\theta}{1+\cos\theta+\epsilon}\rangle$ with
$\epsilon=10^{-5}$; note that when following [26, 12] with $\epsilon=10^{-2}$,
there will be some initial disagreement between network voltage and theory.

The collective dynamics of the network is described by Eqs. (31) with
$I(t)=I_{0}+\kappa\pi R(t)$, which is exact in the thermodynamic limit with
initial conditions $\Phi(0)=\sigma(0)=1$, $\lambda(0)=2$. The initial uniform
voltage distribution amounts to the $\mathcal{M}$-function given by Eq. (38)
with $d=1$. As shown in Fig. 1, the match between simulations (blue) and
theory (orange) is remarkable. Even the initial transients (up to time $t=5$)
towards either of the two bistable attractors is excellently captured by the
theory. By contrast, when restricting the collective dynamics Eqs. (31) to the
two-dimensional OA manifold by setting $\lambda\equiv 0$, with initial
condition $\Phi(0)=\pi R(0)-iV(0)=-v_{0}$, not only are the complex initial
transients no longer captured, but also the dynamics (green curve in Fig. 1)
runs into another attractor different from the actual network dynamics. Note
that in the OA dynamics, the width parameter $d$ of the initial uniform
voltage distribution is lost, so that it is unclear whether the OA dynamics
describes the network evolution from the correct initial condition ($d=1$) or
even from a delta-distribution of voltages ($d=0$).

###  VI.2 Resetting a fraction of neurons induces switching between attractors

As the attentive reader may have noticed, in the example above a slight change
in initial mean voltage, while keeping all the other parameters the same,
results in collective dynamics converging to different attractors: in Fig.
1(a), the dynamics run into a low-activity state; in Fig. 1(b) the dynamics
exhibit a ringing behavior towards a high-activity, yet asynchronous state. As
mentioned above, the attractors of the network dynamics lie on the two-
dimensional OA manifold, allowing for a concise bifurcation analysis of the
bistability between the low- and high-activity states. Below we expand on this
bistable situation, which was already reported in [11] and where the basins of
attraction were explored within the OA manifold. Here we show how our approach
generalizes the basins of attraction in the full state space of QIF spiking
neurons driven by Cauchy white noise and how resetting a fraction of neurons
can induce a switch from one to the other attractor.

As before, we consider global recurrent coupling via instantaneous Dirac
$\delta$-pulses and set $g=0$, so that the total input current is
$I(t)=I_{0}+\kappa\pi R(t)$. In dependence on the scaled parameter
$I_{0}/\Gamma$ one observes a bistability of the stationary firing rate, as
shown in Fig. 2(a). This bistability scenario corresponds to that of Fig. 1(b)
in Ref. [11], where the neurons were assumed heterogeneous ($\Delta=1/4$)
without noise ($\gamma=0$), whereas here we consider identical neurons
($\Delta=0$) with Cauchy noise of strength $\gamma=1/4$. In both cases, we
have the same mean field parameter $\Gamma=\Delta+\gamma=1/4$ and the coupling
strength is $\kappa=\frac{15}{\pi}\sqrt{\Gamma}$.

Figure 2: Attractor switching due to resetting of a fraction of neurons. (a):
Bifurcation analysis yields a region of bistability between two stable,
stationary states (full line) whose basins of attraction on the OA manifold
are separated by the unstable stationary state (dotted line). The pink
vertical line marks the chosen parameter regime $I_{0}/\Gamma=-4$. (b): Basin
boundaries of high-activity attractor: after initializing the system in the
low-activity state, we reset a fraction $\varepsilon$ of oscillators according
to a uniform distribution (blue) or a Cauchy distribution (pink), centered at
$v_{0}$ and having half-width (at half-maximum) 1. Shaded regions indicate
switching to the high-activity state. White domains correspond to decay after
the perturbation back to the initial low-activity state. Model parameters as
in Fig. 1: $\tau_{s}=0$, $g=0$, $\Gamma=0.25$, $I_{0}=-4\Gamma$ and
$\kappa=15\frac{\sqrt{\Gamma}}{\pi}$. Runge-Kutta 4th order method with time
step $dt=10^{-3}$ was used for integration.

We now induce switching from the low- to the high-activity state by resetting
a fraction $\varepsilon$ of the neurons to a predefined voltage distribution.
We consider here identical, noise-driven neurons; this approach cannot be
pursued in the presence of heterogeneity because then the selection of neurons
that are reset crucially influences the collective dynamics. The common input
is fixed at $I_{0}/\Gamma=-4$ in the region of bistability and we initialize
the neurons on the low-activity branch. At the fixed point, the distributions
$W_{0}(v)$ of voltages $v$ and $P_{0}(\theta)$ of phase variables $\theta$ is
Cauchy-Lorentz, or wrapped Cauchy, respectively, as it should be on the OA
manifold. Then we take a portion $0\leq\varepsilon\leq 1$ of the neurons and
redistribute their voltages (phases) according to a new distribution
$W_{1}(v)$ ($P_{1}(\theta)$, correspondingly). Thus, according to expression
(39) we have to start calculations of the transient using full Eqs. (31) with
$\mathcal{M}(k)=(1-\varepsilon)\mathcal{M}_{0}(k)+\varepsilon\mathcal{M}_{1}(k)$.
(Note that in Ref. [36] resettings inside the OA manifold have been
considered).

We use two distributions $W_{1}(v)$ for the resetting:

  1. 1.

A Cauchy-Lorentz distribution (36) with fixed half-width at half-maximum $x=1$
and varying parameter $v_{0}$. Then, $\mathcal{M}_{1}(k)=\frac{\mu k}{1-\mu
k}$ with $\mu=\frac{1-\pi+iv_{0}}{1+\pi-iv_{0}}$.

  2. 2.

A uniform distribution of voltages (37) with half-width $d=1$ and different
$v_{0}$. Because at the initial condition $\sigma(0)=1$, the denominator of
$\mathcal{M}_{1}(-\sigma)$ vanishes, one needs an expansion

$\mathcal{M}_{1}(-\sigma)=-\frac{1+iv_{0}}{2}-\frac{3+\delta^{2}+3v_{0}^{2}}{12}(\sigma-1)+\mathcal{O}((\sigma-1)^{2})$
(46)

to start calculations properly.

Integrating the system of equations (31) with the elaborated constant
functions $\mathcal{M}(k)$, we have observed convergence to one of the
attractors on the OA manifold. The basins of attraction for the different
resetting distributions $W_{1}(v)$ are depicted in the parameter plane
$(v_{0},\varepsilon)$ in Fig. 2(b). Resetting a large fraction of neurons
according to a Cauchy-Lorentz distribution (red), significantly increases the
chances of attractor switching from the low- to the high-activity state
compared to a resetting according to a uniform distribution (blue).

###  VI.3 Phase-dependent resetting from synchronous to asynchronous states

We now add gap junction coupling $g>0$ while keeping the other parameters as
before. For gap junctions, the input current now reads $I(t)=I_{0}+\kappa\pi
R(t)+gV(t)$. Upon increasing $g$, the high-activity state undergoes a Hopf
bifurcation and becomes a stable limit cycle on the OA manifold; the low-
activity state remains a phase locked, fixed point solution, see Fig. 3(a). We
now fix $g=0.35$ and initialize the system in the oscillatory, synchronous
steady state, $Q(t)=Q(t+T)$ with period $T$, that has bifurcated from the
high-activity branch. At different phases of the collective oscillation, we
reset a portion $\varepsilon$ of oscillators with a uniform distribution of
voltages with mean voltage $v_{0}=0$ and half-width $d=1$. We therefore
consider the function
$\mathcal{M}(k)=(1-\varepsilon)\mathcal{M}_{0}(k)+\varepsilon\mathcal{M}_{1}(k)$
as a combination of the Cauchy $\mathcal{M}_{0}(k)=\frac{Qk}{1-Qk}$ with the
expression (38) for the $\mathcal{M}_{1}(k)$. In this special case ($v_{0}=0$,
$d=1$), $\mathcal{M}_{1}(k)$ is given by

$\mathcal{M}_{1}(k)=\frac{-k}{(1+k)^{2}}\Big{[}(1+k)-\pi/2-2\arctan(k)\Big{]}$
(47)

The basin boundary of the asynchronous low-activity state is depicted in blue
in Fig. 3(b), and sensitively depends on the collective oscillation phase and
on the reset fraction $\varepsilon$ of neurons. By contrast, resetting
according to a Cauchy-Lorentz distribution $W_{1}(v)$ centered at 0 with half-
width at half-maximum 1 never induces a switch (for any portion $\varepsilon$
and at any phase of collective oscillation).

Figure 3: Phase-dependent attractor switching from synchronous to asynchronous
states. (a): Bifurcation analysis yields a region of bistability between an
oscillatory, synchronous state (orange), which emerges from a stable high-
activity state via a supercritical Hopf bifurcation, and a stable low-activity
stationary state (full line). The pink vertical line marks the chosen
parameter regime $g=0.35$. (b): Basin boundaries of low-activity attractor:
after initializing the system in synchronous state, we reset a fraction
$\varepsilon$ of oscillators according to a uniform distribution centered at
$0$ and having half-width 1 for every phase of the collective oscillation.
Shaded blue regions correspond to values that induce a switch to the fixed
point. White domain corresponds to decay after the perturbation back to the
collective oscillation state. Model parameters are: $\tau_{s}=0$, $g=0.35$,
$\Gamma=0.25$, $I_{0}=-4\Gamma$ and $\kappa=15\frac{\sqrt{\Gamma}}{\pi}$.
Runge-Kutta 4th order method with time step $dt=10^{-3}$ was used for
integration.

##  VII Discussion & Conclusion

In this paper, we have put forward a finite-dimensional description of large
networks of globally coupled spiking neurons that are described by the
Ermentrout-Kopell canonical model of excitable neuronal systems. Each neuron
can equivalently be expressed as a phase model (theta neuron) or as a
quadratic integrate-and-fire (QIF) neuron with threshold and reset going to
$\pm\infty$, respectively. In the presence of heterogeneous input currents (or
natural frequencies) and/or Cauchy white noise, our formalism is exact in the
thermodynamic limit. The derivation of the set(s) of three complex ordinary
differential equations, (23) or (31), does not rely on assumptions of weak
coupling, separation of time scales, averaging, or any other approximation.
Rather, the assumptions underlying the validity of the low-dimensional
description are that (i) the neurons are all-to-all connected, (ii) noise is
Cauchy (and not Gaussian), and (iii) inputs are distributed according to a
Cauchy-Lorentz distribution; we will comment on these assumptions further
below. We note that in the finite-dimensional reduction, both
situations—neurons are subject to Cauchy white noise, or they receive Cauchy-
Lorentz distributed (time-independent) inputs—result in identical mean-field
dynamics. Yet, only in the former case there is a simple unique correspondence
between the mean-field dynamics via the order parameters and the phase
(voltage) distribution of theta (QIF) neurons. In the case of heterogeneity,
one can calculate the order parameters from the distribution of phases only
under certain analyticity assumptions, e.g., that the density admits an
analytic continuation in the upper complex plane of inputs (or frequencies in
the phase description). Hence, the results of this paper are fully applicable
to noisy ensembles, but some approaches (e.g., the resetting example) are not
suitable for neurons with distributed inputs/heterogeneous frequencies.

We remark that the six-dimensional dynamical reduction presents an important
extension of previous results [8, 9, 10, 11, 12] which were restricted to the
OA (Lorentzian) manifold. According to these results, the collective dynamics
can only be described in the time-asymptotic limit, $t\to\infty$, or if the
initial state of the neurons is meticulously instantiated. By contrast, our
approach allows us to faithfully capture the network dynamics from arbitrary
initial conditions. Furthermore, we can track how the collective dynamics is
eventually attracted to the two-dimensional, attractive OA manifold. We have
derived the global stability of the OA manifold (in the weak sense), whose
attractiveness has already been argued in the literature [33, 34, 35, 37, 38,
15]. In contrast to networks of conventional Kuramoto-type oscillators, where
the convergence rate is given by the degree of heterogeneity and/or noise
strength [15], here we cannot indicate an exact lower bound on the convergence
rate towards the OA manifold; note, however, that here we do not consider
conventional Kuramoto-type oscillators with additive noise and identical
forcing fields $h_{j}\neq h$, but rather parameter-dependent oscillatory
systems with multiplicative noise and oscillator-dependent natural frequencies
$\omega_{j}\neq\omega_{k}$ and forcing fields $h_{j}\neq h_{k}$.

Our theory is valid for any mean-field coupling. As particular examples we
considered all-to-all coupling via instantaneous chemical synapses as well as
via electrical synapses through so-called gap junctions. In Section III we
provide a general framework how chemical interactions via a variety of pulses,
emitted from the pre- to the postsynaptic neuron, can be incorporated through
mean-field variables in our low-dimensional description. While we have focused
here on symmetric pulse profiles, in future work we will investigate the
effect of asymmetric pulses on the collective dynamics of spiking neurons.

As to the specific noisy and quenched inputs, we remark that in case of
Gaussian white noise additional terms appear in the equations for the order
parameters (13), which does not allow for truncating the infinite system (22)
for the $\beta_{n}$. Nevertheless, truncation might yield an approximative
finite-dimensional description similar to that proposed in the Supplementary
Material of [38] and in [39, 40]; this is a subject of a forthcoming research.

Relaxing the nature of quenched heterogeneity appears straightforward, in
particular if the input parameters $\eta$ are drawn from a non-singular
distribution with a finite set of poles outside the real axis, see e.g., Eq.
(2) in [41]. The proposed distribution allows one to approximate both Gaussian
[42] as well as uniform heterogeneity [43, 44] with a finite-dimensional
extension of the OA dynamics to arbitrary accuracy. We leave the corresponding
extension of our exact low-dimensional description of spiking neurons subject
to heterogeneous inputs with $q$-Gaussian or rational distributions,
respectively, for future work. As a side note, we advise caution when dealing
with distributions of inputs that do not comply with the analyticity
assumptions mentioned above. In those cases, transient dynamics can become
nontrivial and the basins of attraction of time-asymptotic solutions depend on
the choice of initial conditions and cannot be treated within the theory
presented here, see, e.g., [45, 46].

Finally, we mention that our approach will hold for more general networks of
spiking neurons. Possible extensions include networks with distributed
synaptic weights $\kappa$, in addition to distributed inputs $\eta$, or
interacting populations of excitatory and inhibitory neurons. Moreover,
synapses can follow more complex synaptic kinetics and can be modeled as
conductances with reversal potentials, that can also be distributed. In these
cases, it may be more important to faithfully capture transient dynamics off
the OA manifold, which can readily be achieved with our exact low-dimensional
description.

###### Acknowledgements.

We thank R. Toenjes and P. Clusella for useful discussions. BP has received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skodowska-Curie grant agreement No 101032806. AP and
RC were supported by the DFG (Grant PI 220/21-1).

## References

  * Strogatz [2003] S. H. Strogatz,  _Sync: The Emerging Science of Spontaneous Order_ (Hyperion, NY, 2003).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization. A Universal Concept in Nonlinear Sciences._ (Cambridge University Press, Cambridge, 2001).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77, 137 (2005).
  * Pikovsky and Rosenblum [2015] A. Pikovsky and M. Rosenblum, Dynamics of globally coupled oscillators: Progress and perspectives, Chaos 25, 097616 (2015).
  * Buzsáki [2006] G. Buzsáki,  _Rhythms of the brain_ (Oxford UP, Oxford, 2006).
  * Izhikevich [2007] E. M. Izhikevich,  _Dynamical Systems in Neuroscience_ (MIT Press, Cambridge, Mass., 2007).
  * Ott and Antonsen [2008] E. Ott and T. M. Antonsen, Low dimensional behavior of large systems of globally coupled oscillators, Chaos 18, 037113 (2008).
  * 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 (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 (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ó  _et al._ [2015] E. Montbrió, D. Pazó, and A. Roxin, Macroscopic description for networks of spiking neurons, Phys. Rev. X 5, 021028 (2015).
  * Laing [2015] C. R. Laing, Exact neural fields incorporating gap junctions, SIAM J. Appl. Dyn. Syst. 14, 1899 (2015).
  * Ashwin  _et al._ [2016] P. Ashwin, S. Coombes, and R. Nicks, Mathematical frameworks for oscillatory network dynamics in neuroscience, J. Math. Neurosci. 6, 1 (2016).
  * Bick  _et al._ [2020] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review, J. Math. Neurosci. 10, 9 (2020).
  * Cestnik and Pikovsky [2022a] R. Cestnik and A. Pikovsky, Exact finite-dimensional reduction for a population of noisy oscillators and its link to Ott-Antonsen and Watanabe-Strogatz theories, arXiv:2207.02302 (2022a).
  * Ermentrout and Kopell [1986] G. B. Ermentrout and N. Kopell, Parabolic bursting in an excitable system coupled with a slow oscillation, SIAM J. Appl. Math. 46, 233 (1986).
  * Ermentrout [1996] B. Ermentrout, Type I membranes, phase resetting curves, and synchrony, Neural Comput. 8, 979 (1996).
  * Kuramoto [1975] Y. Kuramoto,  _in International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki, Vol. 39 (Springer, New York, 1975) p. 420.
  * Sakaguchi and Kuramoto [1986] H. Sakaguchi and Y. Kuramoto, A soluble active rotator model showing phase transition via mutual entrainment, Prog. Theor. Phys. 76, 576 (1986).
  * Cessac and Viéville [2008] B. Cessac and T. Viéville, On dynamics of integrate-and-fire neural networks with conductance based synapses, Front. Comput. Neurosci. , 2 (2008).
  * Kevrekidis  _et al._ [2017] P. G. Kevrekidis, C. I. Siettos, and Y. G. Kevrekidis, To infinity and some glimpses of beyond, Nat. Commun. 8, 1562 (2017).
  * Zygmund [2002] A. Zygmund,  _Trigonometric series_ (Cambridge University Press, 2002).
  * Tönjes and Pikovsky [2020] R. Tönjes and A. Pikovsky, Low-dimensional description for ensembles of identical phase oscillators subject to Cauchy noise, Phys. Rev. E 102, 052315 (2020).
  * Clusella and Montbrió [2022] P. Clusella and E. Montbrió, Regular and sparse neuronal synchronization are described by identical mean field dynamics, arXiv:2208.05515 (2022).
  * Tanaka [2020] T. Tanaka, Low-dimensional dynamics of phase oscillators driven by Cauchy noise, Phys. Rev. E 102, 042220 (2020).
  * Ermentrout [2006] B. Ermentrout, Gap junctions destroy persistent states in excitatory networks, Phys. Rev. E 74, 031918 (2006).
  * Pietras  _et al._ [2019] B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Montbrió, Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks, Phys. Rev. E 100, 042412 (2019).
  * Ariaratnam and Strogatz [2001] J. T. Ariaratnam and S. H. Strogatz, Phase diagram for the Winfree model of coupled nonlinear oscillators, Phys. Rev. Lett. 86, 4278 (2001).
  * Gallego  _et al._ [2017] R. Gallego, E. Montbrió, and D. Pazó, Synchronization scenarios in the Winfree model of coupled oscillators, Phys. Rev. E 96, 042208 (2017).
  * Kato and Jones [2015] S. Kato and M. C. Jones, A tractable and interpretable four-parameter family of unimodal distributions on the circle, Biometrika 102, 181 (2015).
  * Watanabe and Strogatz [1993] S. Watanabe and S. H. Strogatz, Integrability of a globally coupled oscillator array, Phys. Rev. Lett. 70, 2391 (1993).
  * Watanabe and Strogatz [1994] S. Watanabe and S. H. Strogatz, Constants of motion for superconducting Josephson arrays, Physica D 74, 197 (1994).
  * 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 oscillator systems” [Chaos 19, 023117 (2009)], Chaos 21, 025112 (2011).
  * Pietras and Daffertshofer [2016] B. Pietras and A. Daffertshofer, Ott-Antonsen attractiveness for parameter-dependent oscillatory systems, Chaos 26, 103101 (2016).
  * Sarkar and Gupta [2022] M. Sarkar and S. Gupta, Synchronization in the Kuramoto model in presence of stochastic resetting, arXiv:2203.00339 (2022).
  * Engelbrecht and Mirollo [2020] J. R. Engelbrecht and R. Mirollo, Is the Ott-Antonsen manifold attracting?, Phys. Rev. Research 2, 023057 (2020).
  * Cestnik and Pikovsky [2022b] R. Cestnik and A. Pikovsky, Hierarchy of exact low-dimensional reductions for populations of coupled oscillators, Phys. Rev. Lett. 128, 054101 (2022b).
  * Tyulkina  _et al._ [2018] I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, and A. Pikovsky, Dynamics of noisy oscillator populations beyond the Ott-Antonsen ansatz, Phys. Rev. Lett. 120, 264101 (2018).
  * Goldobin  _et al._ [2021] D. S. Goldobin, M. di Volo, and A. Torcini, Reduction methodology for fluctuation driven population dynamics, Phys. Rev. Lett. 127, 038301 (2021).
  * Lafuerza  _et al._ [2010] L. F. Lafuerza, P. Colet, and R. Toral, Nonuniversal results induced by diversity distribution in coupled excitable systems, Phys. Rev. Lett. 105, 084101 (2010).
  * Pyragas and Pyragas [2022] V. Pyragas and K. Pyragas, Mean-field equations for neural populations with $q$-Gaussian heterogeneities, Phys. Rev. E 105, 044402 (2022).
  * Skardal [2018] P. S. Skardal, Low-dimensional dynamics of the Kuramoto model with rational frequency distributions, Phys. Rev. E 98, 022207 (2018).
  * Pietras  _et al._ [2018] B. Pietras, N. Deschle, and A. Daffertshofer, First-order phase transitions in the Kuramoto model with compact bimodal frequency distributions, Phys. Rev. E 98, 062219 (2018).
  * Ott  _et al._ [2008] E. Ott, J. H. Platig, T. M. Antonsen, and M. Girvan, Echo phenomena in large systems of coupled oscillators, Chaos 18, 037115 (2008).
  * Pikovsky and Rosenblum [2011] A. Pikovsky and M. Rosenblum, Dynamics of heterogeneous oscillator ensembles in terms of collective variables, Physica D (Amsterdam) 240, 872 (2011).