# Periodic solutions in next generation neural field models

Carlo R. Laing<EMAIL_ADDRESS>School of Mathematical and
Computational Sciences, Massey University, Private Bag 102-904 NSMC, Auckland,
New Zealand  Oleh E. Omel’chenko<EMAIL_ADDRESS>Institute of
Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25,
14476 Potsdam, Germany.

###### Abstract

We consider a next generation neural field model which describes the dynamics
of a network of theta neurons on a ring. For some parameters the network
supports stable time-periodic solutions. Using the fact that the dynamics at
each spatial location are described by a complex-valued Riccati equation we
derive a self-consistency equation that such periodic solutions must satisfy.
We determine the stability of these solutions, and present numerical results
to illustrate the usefulness of this technique. The generality of this
approach is demonstrated through its application to several other systems
involving delays, two-population architecture and networks of Winfree
oscillators.

neural field model, Riccati equation, theta neuron, Ott/Antonsen, self-
consistency

## Introduction

The collective behaviour of large networks of neurons is a topic of ongoing
interest. One of the simplest forms of behaviour is a periodic oscillation,
which manifests itself as a macroscopic rhythm created by the synchronous
firing of many neurons. Such oscillations have relevance to rhythmic movement
[36], epilepsy [20, 38], schizophrenia [57] neural communication [50] and
EEG/MEG modelling [9], among others. In different networks, oscillations may
arise from mechanisms such as synaptic delays [14], the interaction of
excitatory and inhibitory populations [7, 52], having sufficient connectivity
in an inhibitory network [58], or the finite width of synaptic pulses emitted
by neurons [49]. The modelling and simulation of such networks is essential in
order to investigate their dynamics.

Among the many types of model neurons used when studying networks of neurons,
theta neurons [16], Winfree oscillators [2] and quadratic integrate-and-fire
(QIF) neurons [34] are some of the simplest. These three types of model
neurons have the advantage that their mathematical form often allows infinite
networks of such neurons to be analysed exactly using the Ott-Antonsen method
[44, 45]. We continue along those lines in this paper.

Here we largely consider a spatially-extended network of neurons, in which the
neurons can be thought of as lying on a ring. Such ring networks have been
studied in connection with modelling head direction [61] and working memory
[18, 60], for example, and been studied by others [17, 24, 22]. We consider a
network of theta neurons. The theta neuron is a minimal model for a neuron
which undergoes a saddle-node-on-invariant-circle (SNIC) bifurcation as a
parameter is varied [16]. The theta neuron is exactly equivalent to a
quadratic integrate-and-fire (QIF) neuron, under the assumption of infinite
firing threshold and reset values [37, 4, 15]. The coupling in the network is
nonlocal synaptic coupling, implemented using a spatial convolution with a
translationally-invariant coupling kernel.

We studied such a model in the past [40], concentrating on describing
spatially-uniform states and also stationary “bump” states in which there is a
spatially-localised group of active neurons while the remainder of the network
is quiescent. We determined the existence and stability of such states and
also found regions of parameter space in which neither of these types of
states were stable. Instead, we sometimes found solutions which were periodic
in time. In this paper we study such periodic solutions using a recently-
developed technique [43, 39] which is significantly faster than the standard
approach.

This technique was successfully applied to describe traveling and breathing
chimera states in nonlocally coupled phase oscillators [43, 39]. In this
paper, we generalize this technique to neural models and illustrate its
possibilities with several examples. We start with a ring network of theta
neurons and consider periodic bump states there. We describe a continuation
algorithm, perform linear stability analysis of such states and derive some
useful formulas, for example, for average firing rates. We show that the same
approach works in the presence of delays and for two-population models.
Finally, we show that Winfree oscillators are also treatable by the proposed
technique.

The structure of the paper is as follows. In Sec. I we present the discrete
network model and its continuum-level description. The bulk of the paper is in
Sec II, where we show how to describe periodic states in a self-consistent
way, and show numerical results from implementing our algorithms. Other models
are considered in Sec. III and we end in Sec. IV. The Appendix contains useful
results regarding the complex Riccati equation.

##  I Theta neuron network model

The model we consider first is that in [40, 25, 31], which we briefly present
here. The discrete network consists of $N$ synaptically coupled theta neurons
described by

$\frac{\displaystyle d\theta_{j}}{\displaystyle
dt}=1-\cos\theta_{j}+(1+\cos\theta_{j})(\eta_{j}+\kappa I_{j}),\quad
j=1,\dots,N,$ (1)

where each $\theta_{j}\in[0,2\pi]$ is an angular variable. The constant
$\kappa$ is the overall strength of coupling within the network, and the
current entering the $j$th neuron is $\kappa I_{j}$ where

$I_{j}(t)=\frac{\displaystyle 2\pi}{\displaystyle
N}\sum\limits_{k=1}^{N}K_{jk}P_{n}(\theta_{k}(t))$ (2)

where

$P_{n}(\theta)=a_{n}(1-\cos\theta)^{n}$

is a pulsatile function with a maximum at $\theta=\pi$ (when the neuron fires)
and $a_{n}$ is chosen according to the normalization condition

$\int_{0}^{2\pi}P_{n}(\theta)d\theta=2\pi.$

Increasing $n$ makes $P_{n}(\theta)$ “sharper” and more pulse-like. The
excitability parameters $\eta_{j}$ are chosen from a Lorentzian distribution
with mean $\eta_{0}$ and width $\gamma>0$

$g(\eta)=\frac{\displaystyle\gamma}{\displaystyle\pi}\frac{\displaystyle
1}{\displaystyle(\eta-\eta_{0})^{2}+\gamma^{2}}.$

The connectivity within the network is given by the weights $K_{jk}$ which are
defined by $K_{jk}=K(2\pi(j-k)/N)$ where the coupling kernel is

$K(x)=\frac{\displaystyle 1}{\displaystyle 2\pi}(1+A\cos x)$ (3)

for some constant $A$. Note that the form of coupling implies that the neurons
are equally-spaced around a ring, with periodic boundary conditions. Such a
network can support solutions which are — at a macroscopic level — periodic in
time; see Fig. 3 in [40]. Such solutions are unlikely to be true periodic
solutions of (1), since for a typical realisation of the $\eta_{j}$, one or
more neurons will have extreme values of this parameter, resulting in them not
frequency-locking to other neurons.

Note that the model we study here has only one neuron at each spatial
position, and for $A>0$ the connections between nearby neurons are more
positive than those between distant neurons. For $A>1$ neurons on opposite
sides of the domain inhibit one another, as the connections between them are
negative, giving a “Mexican-hat” connectivity. For $A<0$ connections between
neurons on opposite sides of the domain are more positive than those between
nearby neurons. Such a model with one population of neurons and connections of
mixed sign can be thought of as an approximation of a network with populations
of both excitatory and inhibitory neurons [48, 17].

Using the Ott/Antonsen ansatz [44, 45], one can show that in the limit
$N\to\infty$, the long-term dynamics of the network (1) can be described by

$\displaystyle\frac{\displaystyle\partial z}{\displaystyle\partial t}$
$\displaystyle=$
$\displaystyle\frac{\displaystyle(i\eta_{0}-\gamma)(1+z)^{2}-i(1-z)^{2}}{\displaystyle
2}$ (4) $\displaystyle+$ $\displaystyle\kappa\frac{\displaystyle
i(1+z)^{2}}{\displaystyle 2}\mathcal{K}H_{n}(z),$

where

$(\mathcal{K}\varphi)(x)=\int_{0}^{2\pi}K(x-y)\varphi(y)dy$ (5)

is the convolution of $K$ and $\varphi$ and

$H_{n}(z)=a_{n}\left[C_{0}+\sum\limits_{q=1}^{n}C_{q}\left(z^{q}+\overline{z}^{q}\right)\right],$

where

$C_{q}=\sum\limits_{k=0}^{n}\sum\limits_{m=0}^{k}\frac{\displaystyle\delta_{k-2m,q}(-1)^{k}n!}{\displaystyle
2^{k}(n-k)!m!(k-m)!}.$

For our computations we set $n=2$, so that

$H_{2}(z)=(2/3)[3/2-(z+\bar{z})+(z^{2}+\bar{z}^{2})/4]$

where overline indicates the complex conjugate. Periodic solutions like those
studied below were also found with $n=5$, for example, so the choice of $n$ is
not critical. The wider question of the effects of pulse shape and duration is
an interesting one [47].

Eq. (4) is an integro-differential equation for a complex-valued function
$z(x,t)$, where $x\in[0,2\pi]$ is position on the ring. $z(x,t)$ is a local
order parameter and can be thought of as the average of $e^{i\theta}$ for
neurons in a small neighbourhood of position $x$. The magnitude of $z$ is a
measure of how synchronised the neurons are, whereas its argument gives the
most likely value of $\theta$ [25]. Using the equivalence of theta and QIF
neurons, one can also provide a relevant biological interpretation of $z$.
Namely, defining $W\equiv(1-\overline{z})/(1+\overline{z})$, one can show [27,
37] that $\pi^{-1}\mathrm{Re}\>W$ is the flux through $\theta=\pi$ or the
instantaneous firing rate of neurons at position $x$ and time $t$. Similarly,
if $V_{j}=\tan(\theta_{j}/2)$ is the voltage of the $j$th QIF neuron, the mean
voltage at position $x$ and time $t$ is given by $\mathrm{Im}\>W$. Note that
physically relevant solutions of Eq. (4) must assume values $|z|\leq 1$. In
other words, we are interested only in solutions $z\in\overline{\mathbb{D}}$,
where

$\mathbb{D}=\\{z\in\mathbb{C}\>:\>|z|<1\\}$

is the unit disc in the complex plane.

Equations of the form (4)–(5) are sometimes referred to as next generation
neural field models [8, 11] as they have the form of a neural field model (an
integro-differential equation for a macroscopic quantity such as “activity”
[33, 1]) but are derived exactly from a network like (1), rather than being of
a phenomenological nature.

##  II Periodic states

In this paper we focus on states with periodically oscillating macroscopic
dynamics. For the mean field equation (4) these correspond to periodic
solutions

$z=a(x,t)$

where $a(x,t+T)=a(x,t)$ for some $T>0$. The frequency corresponding to $T$ is
denoted by $\omega=2\pi/T$. An example of a periodic solution is shown in Fig.
1, panels (a) and (b). Fig. 1 (c) shows a realization of the same solution in
the discrete network (1). We note that this solution has a spatio-temporal
symmetry: it is invariant under a time shift of half a period followed by a
reflection about $x=\pi$. However, we do not make use of this symmetry in the
following calculations.

Figure 1: A typical periodic solution of Eq. (4). (a): $\arg{(z(x,t))}$. (b):
$|z(x,t)|$. (c) A realization of this solution in a network of $N=4096$ theta
neurons described by Eq. (1). $\theta_{j}$ is shown in color. Parameters:
$A=-5$, $\eta_{0}=-0.7$, $\kappa=1$, $\gamma=0.01$.

The appearance of periodic solutions in this model is in contrast with
classical one-population neural field models for which they do not seem to
occur [32]. This is another example of next generation neural field models
showing more complex time-dependent behaviour than classical ones [31].

A straight-forward way to study a periodic orbit like that in Fig. 1 would be
discretise Eq. (4) on a spatially-uniform grid and approximate the convolution
using matrix/vector multiplication or otherwise, resulting in a large set of
coupled ordinary differential equations. The periodic solution of these could
then be studied using standard techniques [25], but note that the
computational complexity of this would typically scale as $\sim N^{2}$, where
$N$ is the number of spatial points used in the grid. Instead we propose here
an alternative method based on the ideas from [43, 39], which allows us to
perform the same calculations with only $\sim N$ operations. The main
ingredients of this method are explained in Sections II.1 and II.2. They
include the description of the properties of complex Riccati equation and the
derivation of the self-consistency equation for periodic solutions of Eq. (4).
In Section II.3 we explain how the self-consistency equation can be solved in
the case of coupling function (3). Then in Section II.4 we report some
numerical results obtained with our method. In addition, in Section II.5 we
perform a rigorous linear stability analysis of periodic solutions of Eq. (4)
by considering the spectrum of the corresponding monodromy operator. Finally,
in Section II.6, we show how the mean field equation (4) can be used to
predict the average firing rate distribution in the neural network (1).

###  II.1 Periodic complex Riccati equation and Möbius transformation

By the time rescaling $u(x,t)=z(x,t/\omega)$ we can rewrite Eq. (4) in the
form

$\displaystyle\omega\frac{\displaystyle\partial u}{\displaystyle\partial t}$
$\displaystyle=$
$\displaystyle\frac{\displaystyle(i\eta_{0}-\gamma)(1+u)^{2}-i(1-u)^{2}}{\displaystyle
2}$ (6) $\displaystyle+$ $\displaystyle\kappa\frac{\displaystyle
i(1+u)^{2}}{\displaystyle 2}\mathcal{K}H_{n}(u),$

such that the above $T$-periodic solution of Eq. (4) corresponds to a
$2\pi$-periodic solution of Eq. (6). Then, dividing (6) by $\omega$ and
reordering the terms we can see that Eq. (6) is equivalent to a complex
Riccati equation

$\displaystyle\frac{\displaystyle\partial u}{\displaystyle\partial t}$
$\displaystyle=$ $\displaystyle i\left(W(x,t)+\zeta-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)+2i\left(W(x,t)+\zeta+\frac{\displaystyle
1}{\displaystyle 2\omega}\right)u$ (7) $\displaystyle+$ $\displaystyle
i\left(W(x,t)+\zeta-\frac{\displaystyle 1}{\displaystyle
2\omega}\right)u^{2},$

with the $2\pi$-periodic in $t$ coefficient

$W(x,t)=\frac{\displaystyle\kappa}{\displaystyle 2\omega}\mathcal{K}H_{n}(u)$
(8)

and

$\zeta=\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle 2\omega}.$ (9)

In [43] it was shown (see also Proposition 2 and Remark 2 in the Appendix for
additional details) that independent of the choice of the real-valued periodic
function $W(x,t)$, parameters
$\zeta\in\mathbb{C}_{\mathrm{up}}=\\{z\in\mathbb{C}\>:\>\mathrm{Im}\>z>0\\}$
and $\omega>0$, for every fixed $x\in[0,2\pi]$ Eq. (7) has a unique stable
$2\pi$-periodic solution $U(x,t)$ that lies entirely in the open unit disc
$\mathbb{D}$. Denoting the corresponding solution operator by

$\mathcal{U}\>:\>C_{\mathrm{per}}([0,2\pi];\mathbb{R})\times\mathbb{C}_{\mathrm{up}}\times(0,\infty)\to
C_{\mathrm{per}}([0,2\pi];\mathbb{D}),$

we can write the $2\pi$-periodic solution of interest as

$U(x,t)=\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right).$ (10)

Note that $C_{\mathrm{per}}([0,2\pi];\mathbb{R})$ here denotes the space of
all real-valued continuous $2\pi$-periodic functions, while the notation
$C_{\mathrm{per}}([0,2\pi];\mathbb{D})$ stands for the space of all complex
continuous $2\pi$-periodic functions with values in the open unit disc
$\mathbb{D}$. Importantly, the variable $x$ appears in formula (10) as a
parameter so that the function $W(x,\cdot)\in
C_{\mathrm{per}}([0,2\pi];\mathbb{R})$ with a fixed $x$ is considered as the
first argument of the operator $\mathcal{U}$.

As for the operator $\mathcal{U}$, although it is not explicitly given, its
value can be calculated without resource-demanding iterative methods by
solving exactly four initial value problems for Eq. (7). The rationale for
this approach can be found in [43, Section 4] and is repeated for completeness
in Remark 3 in Appendix. Below we describe its concrete implementation in the
case of formula (10).

We assume that the spatial domain $[0,2\pi]$ is discretised with $N$ points,
$x_{j},j=1,2,\dots N$ and that the functions $U(x,t)$ and $W(x,t)$ are
replaced with their grid counterparts $u_{j}(t)=U(x_{j},t)$ and
$w_{j}(t)=W(x_{j},t)$, respectively. Then, given a set of functions
$w_{j}(t)$, we calculate the corresponding functions $u_{j}(t)$ by performing
the following four steps.

(i) We (somewhat arbitrarily) choose three initial conditions
$u_{j}^{1}(0)=-0.95$, $u_{j}^{2}(0)=0$, $u_{j}^{3}(0)=0.95$ and solve Eq. (7)
with these initial conditions to obtain solutions $u_{j}^{k}(t)$, $j=1,\dots
N$; $k=1,2,3$. In [43] (see also Proposition 2 in the Appendix) it is shown
that for $\gamma>0$ these solutions lie in the open unit disc $\mathbb{D}$.

(ii) Since at each point in space the Poincaré map of Eq. (7) with
$2\pi$-periodic coefficients coincides with a Möbius map $\mathcal{M}(u)$ (see
[43] for detail), we can use the relations
$u_{j}^{k}(2\pi)=\mathcal{M}_{j}(u_{j}^{k}(0))$, $k=1,2,3$, to reconstruct
these maps $\mathcal{M}_{j}$, $j=1,\dots N$. The corresponding formulas are
given in [43, Section 4].

(iii) Now that the Möbius maps $\mathcal{M}_{j}(u)$ are known, their fixed
points can be found by solving $u_{j}^{\ast}=\mathcal{M}_{j}(u_{j}^{\ast})$
for each $j$. This equation is equivalent to a complex quadratic equation,
therefore in general it has two solutions in the complex plane. For
$\gamma>0$, only one of these solutions lies in the unit disc $\mathbb{D}$.

(iv) Using the fixed points $u_{j}^{\ast}\in\mathbb{D}$ of the Möbius maps
$\mathcal{M}_{j}(u)$ as an initial condition in Eq. (7) (i.e. setting
$u(x_{j},0)=u_{j}^{\ast}$) and integrating (7) for a fourth time to $t=2\pi$
we obtain the grid counterpart of a $2\pi$-periodic solution, $U(x,t)$, that
lies entirely in the unit disc $\mathbb{D}$.

We now use this result to show how to derive a self-consistency equation, the
solution of which allows us to determine a $2\pi$-periodic solution of Eq.
(6).

###  II.2 Self-consistency equation

Supposing that Eq. (6) has a $2\pi$-periodic solution, then using formula (8)
we can calculate the corresponding function $W(x,t)$. On the other hand, using
formula (10) we can recover $u(x,t)$. Then the new and the old expressions of
$u(x,t)$ will agree with each other if and only if the function $W(x,t)$
satisfies a self-consistency equation

$W(x,t)=\frac{\displaystyle\kappa}{\displaystyle
2\omega}\mathcal{K}H_{n}\left(\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right)\right),$ (11)

obtained by inserting (10) into (8). In the following we consider Eq. (11) as
a separate equation which must be solved with respect to $W(x,t)$ and
$\omega$. Note that the unknown field $W(x,t)$ has a problem-specific meaning:
It is proportional to the current entering a neuron at position $x$ at time
$t$ due to the activity of all other neurons in the network. The use of self-
consistency arguments to study infinite networks of oscillators goes back to
Kuramoto [23, 56, 55], but such approaches have always focused on steady
states, whereas we consider periodic solutions here.

Note that from a computational point of view, the self-consistency equation
(11) has many advantages. It allows us to reduce the dimensionality of the
problem at least in the case of special coupling kernels with finite number of
Fourier harmonics (see Section II.3). Moreover, its structure is convenient
for parallelization, since the computations of operator $\mathcal{U}$ at
different points $x$ are performed independently. Finally, the main efficiency
is due to the fact that the computation of $\mathcal{U}$ is performed in non-
iterative way.

In the next proposition, we will prove some properties of the solutions of Eq.
(11), which will be used later in Section II.5.

###### Proposition 1

Let the pair $(W(x,t),\omega)$ be a solution of the self-consistency equation
(11) and let $U(x,t)$ be defined by (10). Then

$\left|\exp\left(\int_{0}^{2\pi}M(x,t)dt\right)\right|<1$

where

$\displaystyle M(x,t)$ $\displaystyle=$ $\displaystyle\frac{\displaystyle
i}{\displaystyle\omega}\left[(\kappa\mathcal{K}H_{n}(U)+\eta_{0}+i\gamma)(1+U(x,t))\right.$
(12) $\displaystyle+$ $\displaystyle\left.1-U(x,t)\right].$

Proof: For every fixed $x\in[0,2\pi]$, the function $U(x,t)$ yields a stable
$2\pi$-periodic solution of the complex Riccati equation (7). The
linearization of Eq. (7) around this solution reads

$\frac{\displaystyle dv}{\displaystyle dt}=M(x,t)v,$

where

$\displaystyle M(x,t)$ $\displaystyle=$ $\displaystyle
2i\left(W(x,t)+\zeta+\frac{\displaystyle 1}{\displaystyle 2\omega}\right)$
$\displaystyle+$ $\displaystyle 2i\left(W(x,t)+\zeta-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)U(x,t).$

Moreover, using (8) and (9), we can show that the above expression determines
a function identical to the function $M(x,t)$ in (12). Recalling that $U(x,t)$
is not only a stable but also an asymptotically stable solution of Eq. (7),
see Remark 1 in Appendix, we conclude that the corresponding Floquet
multiplier lies in the open unit disc $\mathbb{D}$. This ends the proof.

###  II.3 Numerical implementation

Eq. (11) describes a periodic orbit, and since Eq. (7) is autonomous we need
to append a pinning condition in order to select a specific solution of Eq.
(11). For a solution of the type shown in Fig. 1 we choose

$\int_{0}^{2\pi}dx\int_{0}^{2\pi}W(x,t)\sin(2t)dt=0.$ (13)

In the following we focus on the case of the cosine coupling (3). It is
straight-forward to show that $(\mathcal{K}H_{n})(x)$ can be written as a
linear combination of $1,\cos{x}$ and $\sin{x}$. However, the system is
translationally invariant in $x$, and we can eliminate this invariance from
Eq. (11) by restricting this equation to its invariant subspace
$\mathrm{Span}\\{1,\sin x\\}$. Then, taking into account that the function
$W(x,t)$ is real, we seek an approximate solution of the system (11), (13)
using a Fourier-Galerkin ansatz

$W(x,t)=\sum\limits_{m=0}^{2F}(v_{m}+w_{m}\sin x)\psi_{m}(t)$ (14)

where $v_{m}$ and $w_{m}$ are real coefficients and $\psi_{m}(t)$ are
trigonometric basis functions

$\displaystyle\psi_{0}(t)$ $\displaystyle=$ $\displaystyle 1,$
$\displaystyle\psi_{m}(t)$ $\displaystyle=$
$\displaystyle\sqrt{2}\cos(nt),\quad\mbox{if}\quad m=2n\quad\mbox{with}\quad
n\in\mathbb{N},$ $\displaystyle\psi_{m}(t)$ $\displaystyle=$
$\displaystyle\sqrt{2}\sin(nt),\quad\mbox{if}\quad m=2n-1\quad\mbox{with}\quad
n\in\mathbb{N}.$

Our typical choice of the number of harmonics in (14) is $F=10$. To exactly
represent $W(x,t)$ in (14) would require an infinite number of terms in the
series, so using a finite value of $F$ introduces an approximation in our
calculations. However, the excellent agreement between our calculations with
$F=10$ and those from full simulations of (4) (shown below) indicate that such
an approximation is justified.

Using the scalar product

$\langle u,v\rangle=\frac{\displaystyle
1}{\displaystyle(2\pi)^{2}}\int_{0}^{2\pi}dx\int_{0}^{2\pi}u(x,t)v(x,t)dt$

we project Eq. (11) on different spatio-temporal Fourier modes to obtain the
system

$\displaystyle v_{m}$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\kappa}{\displaystyle 2\omega}\left\langle
H_{n}\left(\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right)\right),\psi_{m}(t)\right\rangle,$ (15) $\displaystyle
w_{m}$ $\displaystyle=$ $\displaystyle\frac{\displaystyle\kappa
A}{\displaystyle 2\omega}\left\langle
H_{n}\left(\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right)\right),\right.$ (16)
$\displaystyle\left.\vphantom{H_{n}\left(\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right)\right)}\psi_{m}(t)\sin x\right\rangle,$

for $m=0,1,\dots,2F$. Eqs. (15) and (16), together with (13), are a set of
$2(2F+1)+1=4F+3$ equations for the $4F+3$ unknowns
$v_{0},v_{1},\dots,v_{2F},w_{0},w_{1},\dots,w_{2F},\omega$, which must be
solved simultaneously. We solve them using Newton’s method and find
convergence within 3 or 4 iterations.

In simple terms, suppose we have somewhat accurate estimates of
$v_{0},v_{1},\dots,v_{2F},w_{0},w_{1},\dots,w_{2F},\omega$. These can be
inserted into (14) to calculate the function $W(x,t)$. Then one can calculate
a periodic solution of Eq. (7) with the specified $W(x,t)$ by formula (10) and
finally calculate $H_{n}(\mathcal{U})$ and insert this into (15) and (16) and
perform the projections. We want the difference between the new values of
$v_{0},v_{1},\dots,v_{2F},w_{0},w_{1},\dots,w_{2F}$ and the initial values of
them to be zero, and for (13) to hold. This determines the $4F+3$ equations we
need to solve. The solutions of these equations can be followed as a parameter
is varied in the standard way [26]. Note that these calculations involve
discretising the spatial domain with $N$ points. However, the number of
unknowns ($4F+3$) is significantly less than $N$.

###  II.4 Results

The results of following the solution shown in Fig. 1 as $\eta_{0}$ is varied
are shown in Fig. 2, along with values measured from direct simulations of Eq.
(4). The period seems to become arbitrarily large as $\eta_{0}$ approaches
$-2.32$, and the solution approaches a heteroclinic connection, spending more
and more time near two symmetric states which are mapped to one another under
the transformation $x\to-x$. The solution becomes unstable at
$\eta_{0}\approx-0.5$ through a subcritical torus (or secondary Hopf)
bifurcation. (This was determined by finding the Floquet multipliers of the
periodic solution in the conventional way; results not shown.) For these
calculations we used $N=256$ spatial points.

Figure 2: Period of the type of solution shown in Fig. 1 as a function of
$\eta_{0}$ (solid curve). The circles show values measured from direct
simulations of Eq. (4). Other parameters: $A=-5$, $\kappa=1$, $\gamma=0.01$.

For more negative values of $\eta_{0}$ than those shown in Fig. 2, another
type of periodic solution is stable: see Fig. 3. Such a solution does not have
the spatio-temporal symmetry of the solution shown in Fig. 1. However, we can
follow it in just the same way as the parameter $\eta_{0}$ is varied, and we
obtain the results shown in Fig. 4. This periodic orbit appears to be
destroyed in a supercritical Hopf bifurcation as $\eta_{0}$ is decreased
through approximately $-3.2$, and become unstable to a wandering pattern at
$\eta_{0}$ is increased through approximately $-2.34$.

Figure 3: Another periodic solution of Eq. (4). (a): $\arg{(z(x,t))}$. (b):
$|z(x,t)|$. Parameters: $A=-5$, $\eta_{0}=-2.5$, $\kappa=1$, $\gamma=0.01$.
Figure 4: Period of the type of solution shown in Fig. 3 as a function of
$\eta_{0}$ (solid curve). The circles show values measured from direct
simulations of Eq. (4). Other parameters: $A=-5$, $\kappa=1$, $\gamma=0.01$.

Note that the left asymptote in Fig. 2 coincides with the right asymptote in
Fig. 4. On the other hand, we note that two patterns shown in Figs. 1 and 3
have different spatiotemporal symmetries, therefore due to topological reasons
they cannot continuously transform into each other. Similar bifurcation
diagrams where parameter ranges of two patterns with different symmetries are
separated by heteroclinic or homoclinic bifurcations were found for non-
locally coupled Kuramoto-type phase oscillators [42] and seem to be a general
mechanism which, however, needs additional investigation.

###  II.5 Stability of breathing bumps

Given a $T$-periodic solution $a(x,t)$ of Eq. (4), we can perform its linear
stability analysis, using the approach proposed in [39]. Before doing this, we
write

$H_{n}(z)=a_{n}C_{0}+2\mbox{Re}[D_{n}(z)]$

where

$D_{n}(z)=a_{n}\sum_{q=1}^{n}C_{q}z^{q},$

to emphasise that $H_{n}(z)$ is always real. Now, we insert the ansatz
$z(x,t)=a(x,t)+v(x,t)$ into Eq. (4) and linearize the resulting equation with
respect to small perturbations $v(x,t)$. Thus, we obtain a linear integro-
differential equation

$\frac{\displaystyle\partial v}{\displaystyle\partial
t}=\mu(x,t)v+\kappa\frac{\displaystyle i(1+a(x,t))^{2}}{\displaystyle
2}\mathcal{K}\left(D_{n}^{\prime}(a)v+\overline{D_{n}^{\prime}(a)}\overline{v}\right),$
(17)

where

$\mu(x,t)=[i(\eta_{0}+\kappa\mathcal{K}H_{n}(a))-\gamma](1+a(x,t))+i(1-a(x,t))$
(18)

and

$D_{n}^{\prime}(z)=\frac{\displaystyle d}{\displaystyle
dz}D_{n}(z)=a_{n}\sum\limits_{q=1}^{n}qC_{q}z^{q-1}.$

Note that Eq. (17) coincides with the Eq. (5.1) from [40], except that the
coefficients $a(x,t)$ and $\mu(x,t)$ are now time-dependent. Since Eq. (17)
contains the complex-conjugated term $\overline{v}$, it is convenient to
consider this equation along with its complex-conjugate

$\frac{\displaystyle\partial\overline{v}}{\displaystyle\partial
t}=\overline{\mu}(x,t)\overline{v}-\kappa\frac{\displaystyle
i(1+\overline{a}(x,t))^{2}}{\displaystyle
2}\mathcal{K}\left(D_{n}^{\prime}(a)v+\overline{D_{n}^{\prime}(a)}\overline{v}\right).$

This pair of equations can be written in the operator form

$\frac{\displaystyle dV}{\displaystyle dt}=\mathcal{A}(t)V+\mathcal{B}(t)V,$
(19)

where $V(t)=(v_{1}(t),v_{2}(t))^{\mathrm{T}}$ is a function with values in
$C_{\mathrm{per}}([0,2\pi];\mathbb{C}^{2})$, and

$\mathcal{A}(t)V=\left(\begin{array}[]{cc}\mu(x,t)&0\\\\[5.69054pt]
0&\overline{\mu}(x,t)\end{array}\right)\left(\begin{array}[]{c}v_{1}\\\\[5.69054pt]
v_{2}\end{array}\right),$

and

$\mathcal{B}(t)V=\frac{\displaystyle i\kappa}{\displaystyle
2}\left(\begin{array}[]{cc}(1+a(x,t))^{2}&0\\\\[5.69054pt]
0&-(1+\overline{a}(x,t))^{2}\end{array}\right)$
$\times\left(\begin{array}[]{c}\mathcal{K}\left(D_{n}^{\prime}(a)v_{1}\right)\\\\[5.69054pt]
\mathcal{K}\left(\overline{D_{n}^{\prime}(a)}v_{2}\right)\end{array}\right).$
(20)

For every fixed $t$ the operators $\mathcal{A}(t)$ and $\mathcal{B}(t)$ are
linear operators from $C_{\mathrm{per}}([0,2\pi];\mathbb{C}^{2})$ into itself.
Moreover, they both depend continuously on $t$ and thus their norms are
uniformly bounded for all $t\in[0,T]$.

Recall that the question of linear stability of $a(x,t)$ in Eq. (4) is
equivalent to the question of linear stability of the zero solution in Eq.
(17), and hence to the question of linear stability of the zero solution in
Eq. (19). Moreover, using the general theory of periodic differential
equations in Banach spaces, see [13, Chapter V], the last question can be
reduced to the analysis of the spectrum of the monodromy operator
$\mathcal{E}(T)$ defined by the operator exponent

$\mathcal{E}(t)=\exp\left(\int_{0}^{t}(\mathcal{A}(t^{\prime})+\mathcal{B}(t^{\prime}))dt^{\prime}\right).$

The analysis of Eq. (19) in the case when $\mathcal{A}(t)$ is a matrix
multiplication operator and $\mathcal{B}(t)$ is an integral operator similar
to (20) has been performed in [39, Section 4]. Repeating the same arguments we
can demonstrate that the spectrum of the monodromy operator $\mathcal{E}(T)$
is bounded and symmetric with respect to the real axis of the complex plane.
Moreover, it consists of two qualitatively different parts:

(i) the essential spectrum, which is given by the formula

$\sigma_{\mathrm{ess}}=\left\\{\exp\left(\int_{0}^{T}\mu(x,t)dt\right)\>:\>x\in[0,2\pi]\right\\}\cup\\{\mathrm{c.c.}\\}$
(21)

(ii) the discrete spectrum $\sigma_{\mathrm{disc}}$ that consists of finitely
many isolated eigenvalues $\lambda$, which can be found using a characteristic
integral equation, as explained in [39, Section 4].

Note that if $a(x,t)$ is obtained by solving the self-consistency equation
(11) and hence it satisfies

$a(x,t/\omega)=U(x,t)=\mathcal{U}\left(W(x,t),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right),$

where $(W(x,t),\omega)$ is a solution of Eq. (11), then we can use Proposition
1 and formula (18) to show

$\left|\exp\left(\int_{0}^{T}\mu(x,t)dt\right)\right|<1\quad\mbox{for
all}\quad x\in[0,2\pi].$

In this case, the essential spectrum $\sigma_{\mathrm{ess}}$ lies in the open
unit disc $\mathbb{D}$ and therefore it cannot contribute to any linear
instability of the zero solution of Eq. (19).

To illustrate the usefulness of formula (21), in Fig. 5 (a) we plot the
essential spectrum for the periodic solution shown in Fig. 1. In Fig. 5 (b) we
show the Floquet multipliers of the same periodic solution, where we have
found the solution and its stability in the conventional way, of discretizing
the domain and finding a periodic solution of a large set of coupled ordinary
differential equations. In panel (b) we see several real Floquet multipliers
that do not appear in panel (a); these are presumably part of the discrete
spectrum. Note that calculating the discrete spectrum by the method of [39,
Section 4] is numerically difficult, so we do not do that here.

Figure 5: (a) The essential spectrum given by (21) for the periodic solution
shown in Fig. 1. (b) Floquet multipliers of the same periodic solution found
using the technique explained at the start of Sec. II. For both calculations
the spatial domain has been discretized using 512 evenly spaced points.

###  II.6 Formula for firing rates

One quantity of interest in a network of model neurons such as (1) is their
firing rate. The firing rate of the $k$th neuron is defined by

$f_{k}=\frac{\displaystyle 1}{\displaystyle
2\pi}\left\langle\frac{\displaystyle d\theta_{k}}{\displaystyle
dt}\right\rangle_{\mathrm{T}},$

where the angled brackets $\langle\cdot\rangle_{\mathrm{T}}$ indicate a long-
time average. In the case of large $N$, we can also consider the average
firing rate

$f(x)=\frac{\displaystyle
1}{\displaystyle\\#\\{k:|x_{k}-x|<\pi/\sqrt{N}\\}}\sum\limits_{|x_{k}-x|<\pi/\sqrt{N}}f_{k},$
(22)

where $x_{k}=2\pi k/N$ is the spatial positions of the $k$th neuron and the
averaging takes place over all neurons in the $(\pi/\sqrt{N})$-vicinity of the
point $x\in[0,2\pi]$. Note that while the individual firing rates $f_{k}$ are
usually randomly distributed due to the randomness of the excitability
parameters $\eta_{k}$, the average firing rate $f(x)$ converges to a
continuous (and even smooth) function for $N\to\infty$. Moreover, the exact
prediction of the limit function $f(x)$ can be given, using only the
corresponding solution $z(x,t)$ of Eq. (4). To show this, we write Eq. (1) as

$\frac{\displaystyle d\theta_{k}}{\displaystyle
dt}=\mathrm{Re}\left\\{1-e^{i\theta_{k}}+(\eta_{k}+\kappa
I_{k})(1+e^{i\theta_{k}})\right\\}.$ (23)

We recall that in deriving (4) from (1) we introduce a probability
distribution $\rho(\theta,x,\eta,t)$ which satisfies a continuity equation
[27, 41, 25]. At a given time $t$, $\rho(\theta,x,\eta,t)d\theta d\eta dx$ is
the probability that a neuron with a position in $[x,x+dx]$ and intrinsic
drive in $[\eta,\eta+d\eta]$ has its phase in $[\theta,\theta+d\theta]$.
Moreover, in the case of the Lorentzian distribution of parameters $\eta_{k}$,
the probability distribution $\rho(\theta,x,\eta,t)$ satisfies the relations

$\displaystyle\int_{-\infty}^{\infty}d\eta\int_{0}^{2\pi}\rho(\theta,x,\eta,t)e^{i\theta}d\theta=z(x,t),$
(24)
$\displaystyle\int_{-\infty}^{\infty}d\eta\int_{0}^{2\pi}\eta\rho(\theta,x,\eta,t)e^{i\theta}d\theta=(\eta_{0}+i\gamma)z(x,t),$
(25)

which are obtained by a standard contour integration in the complex plane
[44]. The relation (24) has been already used to calculate the continuum limit
analog of (2)

$I(x,t)=\int_{0}^{2\pi}K(x-y)H_{n}(z(y,t))dy=\mathcal{K}H_{n}(z).$

Inserting this instead of $I_{k}(t)$ into Eq. (23) and replacing the time and
index averaging in (22) with the corresponding averaging over the probability
denisty, we obtain

$\displaystyle
f(x)=\lim_{\tau\to\infty}\frac{1}{2\pi\tau}\int_{0}^{\tau}\int_{-\infty}^{\infty}\int_{0}^{2\pi}\mathrm{Re}\left\\{1-e^{i\theta}\right.$
$\displaystyle\left.+(\eta+\kappa
I(x,t))(1+e^{i\theta})\right\\}\rho(\theta,x,\eta,t)d\theta\ d\eta\ dt.$

The two inner integrals in the above formula can be simplified using the
relations (24), (25) and the standard normalization condition for
$\rho(\theta,x,\eta,t)$. Thus we obtain

$\displaystyle
f(x)=\lim_{\tau\to\infty}\frac{1}{2\pi\tau}\int_{0}^{\tau}\mathrm{Re}\left\\{1-z(x,t)\right.$
$\displaystyle\left.+(\eta_{0}+i\gamma+\kappa I(x,t))(1+z(x,t))\right\\}dt.$

Moreover, if $z=a(x,t)$ is a $T$-periodic solution of Eq. (4), then the long-
time average is the same as an average over one period. So, in the periodic
case, the continuum limit average firing rate equals

$\displaystyle f(x)=\frac{\displaystyle 1}{\displaystyle 2\pi
T}\int_{0}^{T}\mathrm{Re}\left\\{1-a(x,t)\right.$
$\displaystyle\left.+(\eta_{0}+i\gamma+\kappa I(x,t))(1+a(x,t))\right\\}dt.$
(26)

(Note that with simple time rescaling formula (26) can be rewritten in terms
of a $2\pi$-periodic solution of the complex Riccati equation (7),
$u(x,t)=a(x,t/\omega)$.) The expression (26) is different from the firing rate
expression given in Sec. I, but both are equally valid.

Results for a pattern like that shown in Fig. 1 are given in Fig. 6, where we
show both $f(x)$ (from (26)) and values extracted from a long simulation of
Eq. (1). The agreement is very good.

Figure 6: Average firing rate for a pattern like that shown in Fig. 1. The
curve shows $f(x)$ as given by (26). The dots show values measured from direct
simulations of Eq. (1). For the discrete simulation, $N=2^{14}$ neurons were
used and the average frequency profile, $\\{f_{j}\\},j=1,2\dots N$, was
convolved with a spatial Gaussian filter of standard deviation $0.01$ before
plotting. For clarity, not all points are shown. Other parameters: $A=-5$,
$\eta_{0}=-0.9$, $\kappa=1$, $\gamma=0.01$.

##  III Other models

We now demonstrate how the approach presented above can be applied to various
other neural models.

###  III.1 Delays

Delays in neural systems are ubiquitous due to the finite velocity at which
action potentials propagate as well as to both dendritic and synaptic
processing [51, 12, 3, 14]. Here we assume that all $I_{j}(t)$ are delayed by
a fixed amount $\tau$, i.e. we have Eq. (1) but we replace (2) by

$I_{j}(t)=\frac{\displaystyle 2\pi}{\displaystyle
N}\sum\limits_{k=1}^{N}K_{jk}P_{n}(\theta_{k}(t-\tau)).$ (27)

The mean field equation is now

$\displaystyle\frac{\displaystyle\partial z}{\displaystyle\partial t}$
$\displaystyle=$
$\displaystyle\frac{\displaystyle(i\eta_{0}-\gamma)(1+z)^{2}-i(1-z)^{2}}{\displaystyle
2}$ (28) $\displaystyle+$ $\displaystyle\kappa\frac{\displaystyle
i(1+z)^{2}}{\displaystyle 2}\mathcal{K}H_{n}(z(x,t-\tau)).$

We can write this equation in the same form as Eq. (7) but now

$W(x,t)=\frac{\displaystyle\kappa}{\displaystyle
2\omega}\mathcal{K}H_{n}(u(x,t-\omega\tau)),$ (29)

and the corresponding self-consistency equation is also time-delayed

$W(x,t)=\frac{\displaystyle\kappa}{\displaystyle
2\omega}\mathcal{K}H_{n}\left(\mathcal{U}\left(W(x,t-\omega\tau),\frac{\displaystyle\eta_{0}+i\gamma}{\displaystyle
2\omega},\omega\right)\right).$

We expand $W(x,t)$ as in (14), and given $W(x,t)$, we find the relevant
$2\pi$-periodic solution of Eq. (7) as above. The only difference comes in
evaluating the projections (15) and (16). Instead of using $H_{n}(U(x,t))$ in
the scalar products, we need to use $H_{n}(U(x,t-\omega\tau))$.

Since $U(x,t)$ is $2\pi$-periodic in time, we can evaluate it at any time
using just its values for $t\in[0,2\pi]$. Specifically,

$U(x,t-\omega\tau)=\begin{cases}U(x,2\pi+t-\omega\tau),&0\leq
t\leq\omega\tau,\\\\[2.84526pt] U(x,t-\omega\tau),&\omega\tau<t\leq
2\pi.\end{cases}$ (30)

Note that this approach would also be applicable if one had a distribution of
delays [35, 30] or even state-dependent delays [21].

As an example, we show in Fig. 7 the results of varying the delay $\tau$ on a
solution of the form shown in Fig. 3. Increasing $\tau$ leads to the
destruction of the periodic solution in an apparent supercritical Hopf
bifurcation.

Figure 7: The vertical axis relates to averaging $W(x,t)$ over $x$. For a
periodic solution, the maximum and minimum over one period of this quantity is
plotted. The black horizontal line corresponds to the steady state which is
stable at $\tau=2.5$. Circles: measured from direct simulations of Eq. (28).
Other parameters: $A=-5$, $\eta_{0}=-2$, $\kappa=1$, $\gamma=0.1$.

###  III.2 Two populations

Neurons fall into two major categories: excitatory and inhibitory. A model
consisting of a single population with a coupling function of the form (3) is
often used as an approximation to a two-population model [17]. Here we
consider a two population model which supports a travelling wave. The mean
field equations are

$\displaystyle\frac{\displaystyle\partial u}{\displaystyle\partial t}$
$\displaystyle=$
$\displaystyle\frac{\displaystyle(i\eta_{u}-\gamma)(1+u)^{2}-i(1-u)^{2}}{\displaystyle
2}$ (31) $\displaystyle+$ $\displaystyle\frac{\displaystyle
i(1+u)^{2}}{\displaystyle
2}\left[w_{\mathrm{ee}}\mathcal{K}H_{n}(u)-w_{\mathrm{ei}}\mathcal{K}H_{n}(v)\right],$
$\displaystyle\frac{\displaystyle\partial v}{\displaystyle\partial t}$
$\displaystyle=$
$\displaystyle\frac{\displaystyle(i\eta_{v}-\gamma)(1+v)^{2}-i(1-v)^{2}}{\displaystyle
2}$ (32) $\displaystyle+$ $\displaystyle\frac{\displaystyle
i(1+v)^{2}}{\displaystyle
2}\left[w_{\mathrm{ie}}\mathcal{K}H_{n}(u)-w_{\mathrm{ii}}\mathcal{K}H_{n}(v)\right]$

where $u(x,t)$ is the complex-valued order parameter for the excitatory
population and $v(x,t)$ is that for the inhibitory population. The non-
negative connectivity kernel between and within populations is the same:

$K(x)=\frac{1}{2\pi}(1+\cos{x})$

and there are four connection strengths within and between populations:
$w_{\mathrm{ee}}$, $w_{\mathrm{ei}}$, $w_{\mathrm{ie}}$ and $w_{\mathrm{ii}}$.
Similar models have been studied in [6, 48].

For some parameter values, such a system supports a travelling wave with a
constant profile. Such a wave can be found very efficiently using the
techniques discussed here, and that was done for a travelling chimera in [43].
However, here we consider a slightly different case: that where the mean drive
to the excitatory population, $\eta_{u}$, is spatially modulated. We thus
write

$\eta_{u}=\eta_{0}+\epsilon\sin{x}.$

For small $|\epsilon|$ the travelling wave persists, but not with a constant
profile. An example is shown in Fig. 8. Note that such a solution is periodic
in time.

Figure 8: A modulated travelling wave solution of Eqs. (31)–(32). (a): $|u|$;
(b): $|v|$. Other parameters: $\eta_{0}=0.1$, $\eta_{v}=0.1$, $\epsilon=0.01$,
$\gamma=0.03$, $w_{\mathrm{ee}}=1$, $w_{\mathrm{ei}}=0.7$,
$w_{\mathrm{ie}}=0.3$, $w_{\mathrm{ii}}=0.1$.

By rescaling time we can write Eqs. (31)–(32) as

$\displaystyle\frac{\displaystyle\partial\tilde{u}}{\displaystyle\partial t}$
$\displaystyle=$ $\displaystyle
i\left(w_{\mathrm{ee}}W_{u}-w_{\mathrm{ei}}W_{v}+\zeta_{u}-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)$ (33) $\displaystyle+$ $\displaystyle
2i\left(w_{\mathrm{ee}}W_{u}-w_{\mathrm{ei}}W_{v}+\zeta_{u}+\frac{\displaystyle
1}{\displaystyle 2\omega}\right)\tilde{u}$ $\displaystyle+$ $\displaystyle
i\left(w_{\mathrm{ee}}W_{u}-w_{\mathrm{ei}}W_{v}+\zeta_{u}-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)\tilde{u}^{2},$
$\displaystyle\frac{\displaystyle\partial\tilde{v}}{\displaystyle\partial t}$
$\displaystyle=$ $\displaystyle
i\left(w_{\mathrm{ie}}W_{u}-w_{\mathrm{ii}}W_{v}+\zeta_{v}-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)$ (34) $\displaystyle+$ $\displaystyle
2i\left(w_{\mathrm{ie}}W_{u}-w_{\mathrm{ii}}W_{v}+\zeta_{v}+\frac{\displaystyle
1}{\displaystyle 2\omega}\right)\tilde{v}$ $\displaystyle+$ $\displaystyle
i\left(w_{\mathrm{ie}}W_{u}-w_{\mathrm{ii}}W_{v}+\zeta_{v}-\frac{\displaystyle
1}{\displaystyle 2\omega}\right)\tilde{v}^{2},$

where $\tilde{u}(x,t)\equiv u(x,t/\omega)$, $\tilde{v}(x,t)\equiv
v(x,t/\omega)$,

$\displaystyle W_{u}(x,t)$ $\displaystyle=$ $\displaystyle\frac{\displaystyle
1}{\displaystyle 2\omega}\mathcal{K}H_{n}(u),$ (35) $\displaystyle W_{v}(x,t)$
$\displaystyle=$ $\displaystyle\frac{\displaystyle 1}{\displaystyle
2\omega}\mathcal{K}H_{n}(v)$ (36)

and

$\zeta_{u}=\frac{\displaystyle\eta_{u}+i\gamma}{\displaystyle
2\omega}=\frac{\displaystyle\eta_{0}+\epsilon\sin{x}+i\gamma}{\displaystyle
2\omega},\qquad\zeta_{v}=\frac{\displaystyle\eta_{v}+i\gamma}{\displaystyle
2\omega}.$

In the same way as above, we can derive self-consistency equations of the form
(11) for $W_{u}(x,t)$ and $W_{v}(x,t)$. Doing so, we obtain a system of two
coupled equations

$\displaystyle W_{u}(x,t)$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\kappa}{\displaystyle
2\omega}\mathcal{K}H_{n}\bigg{(}\mathcal{U}\bigg{(}w_{\mathrm{ee}}W_{u}(x,t)-w_{\mathrm{ei}}W_{v}(x,t),$
$\displaystyle\left.\left.\frac{\displaystyle\eta_{0}+\epsilon\sin{x}+i\gamma}{\displaystyle
2\omega},\omega\right)\right),$ $\displaystyle W_{u}(x,t)$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\kappa}{\displaystyle
2\omega}\mathcal{K}H_{n}\bigg{(}\mathcal{U}\bigg{(}w_{\mathrm{ie}}W_{u}(x,t)-w_{\mathrm{ii}}W_{v}(x,t),$
$\displaystyle\left.\left.\frac{\displaystyle\eta_{v}+i\gamma}{\displaystyle
2\omega},\omega\right)\right).$

One difference between this model and the ones studied above is that the
solution cannot be shifted by a constant amount in $x$ to ensure that it is
always even (or odd) about a particular point in the domain. Thus we need to
write

$\displaystyle W_{u}(x,t)$ $\displaystyle=$
$\displaystyle\sum\limits_{m=0}^{2F}(v_{m}^{u}+w_{m}^{u}\sin x+z_{m}^{u}\cos
x)\psi_{m}(t)$ (37) $\displaystyle W_{v}(x,t)$ $\displaystyle=$
$\displaystyle\sum\limits_{m=0}^{2F}(v_{m}^{v}+w_{m}^{v}\sin x+z_{m}^{v}\cos
x)\psi_{m}(t)$ (38)

These equations contain $6(2F+1)$ unknowns and we find them (and $\omega$) in
the same way as above, by projecting the self-consistency equations for
$W_{u}(x,t)$ and $W_{v}(x,t)$ onto the different spatio-temporal Fourier modes
to obtain equations similar to (15)–(16).

The results of varying the heterogeneity strength $\epsilon$ are shown in Fig.
9. Increasing heterogeneity decreases the period of oscillation, and
eventually the travelling wave appears to be destroyed in a saddle-node
bifurcation.

Figure 9: Period, $T$, of a modulated travelling wave solution of Eqs.
(31)–(32) as a function of heterogeneity strength $\epsilon$. Circles are from
direct simulation of Eqs. (31)–(32). Other parameters are as in Fig. 8.

We conclude this section by noting that for some parameter values the model
(31)–(32) can show periodic solutions which do not travel, like those shown in
Sec. II. Likewise, the model in Sec. I can support travelling waves for
$\kappa=2$.

###  III.3 Winfree oscillators

One of the first models of interacting oscillators studied is the Winfree
model [2, 29, 46, 19]. We consider a spatially-extended network of Winfree
oscillators whose dynamics are given by

$\frac{\displaystyle d\theta_{j}}{\displaystyle
dt}=\omega_{j}+\epsilon\frac{2\pi
Q(\theta_{j})}{N}\sum_{k=1}^{N}K_{jk}P(\theta_{k})$

where $K_{jk}=K(2\pi|j-k|/N)$ for some $2\pi$-periodic coupling function $K$,
$Q(\theta)=-\sin{\theta}/\sqrt{\pi}$ and $P(\theta)=(2/3)(1+\cos{\theta})^{2}$
is a pulsatile function with its peak at $\theta=0$. The $\omega_{j}$ are
randomly chosen from a Lorentzian with centre $\omega_{0}$ and width $\Delta$
and $\epsilon$ is the overall coupling strength.

In the limit $N\to\infty$, using the Ott/Antonsen ansatz, one finds that the
network is described by the equation [28]

$\frac{\displaystyle\partial z}{\displaystyle\partial
t}=\frac{\epsilon}{2\sqrt{\pi}}\mathcal{K}\widehat{H}(z)+(i\omega_{0}-\Delta)z-\frac{\epsilon}{2\sqrt{\pi}}z^{2}\mathcal{K}\widehat{H}(z),$
(39)

where the integral operator $\mathcal{K}$ is again defined by (5) and

$\widehat{H}(z)=(2/3)[3/2+z+\bar{z}+(z^{2}+\bar{z}^{2})/4].$

A typical periodic solution of Eq. (39) for the choice

$K(x)=0.1+0.3\cos{x},$

is shown in Fig. 10.

If we rescale time by the frequency of periodic solution $\omega>0$, defining
$u(x,t)=z(x,t/\omega)$, and denote

$W(x,t)=\frac{\epsilon}{2\sqrt{\pi}\omega}\mathcal{K}\widehat{H}(u),$

then Eq. (39) can be recast as a complex Riccati equation

$\frac{\displaystyle\partial u}{\displaystyle\partial
t}=W(x,t)+i\frac{\displaystyle\omega_{0}+i\Delta}{\displaystyle\omega}u-W(x,t)u^{2}.$
(40)

From the Proposition 2 in [43], it follows that for every $\omega,\Delta>0$,
$\omega_{0}\in\mathbb{R}$ and for every real-valued, $2\pi$-periodic in $t$
function $W(x,t)$, Eq. (40) has a $2\pi$-periodic solution lying in the open
unit disc $\mathbb{D}$. Denoting the corresponding solution operator

$\widehat{\mathcal{U}}\>:\>C_{\mathrm{per}}([0,2\pi];\mathbb{R})\times\mathbb{C}_{\mathrm{up}}\to
C_{\mathrm{per}}([0,2\pi];\mathbb{D}),$

we easily obtain a self-consistency equation for periodic solutions of Eq.
(40)

$W(x,t)=\frac{\epsilon}{2\sqrt{\pi}\omega}\mathcal{K}\widehat{H}\left(\widehat{\mathcal{U}}\left(W(x,t),\frac{\displaystyle\omega_{0}+i\Delta}{\displaystyle\omega}\right)\right).$
(41)

Since the Poincaré map of Eq. (40) coincides with the Möbius transformation,
we can again use the calculation scheme of Sec. II.1 to find the value of
operator $\widehat{\mathcal{U}}$. Thus we can solve Eq. (41) numerically for
the real-valued field $W(x,t)$ and frequency $\omega$.

Figure 10: A periodic solution of Eq. (39). (a): $\arg{(z(x,t))}$. (b):
$|z(x,t)|$. Parameters: $\omega_{0}=1$, $\Delta=0.1$, $\epsilon=2.1$.

Following the periodic solution shown in Fig. 10 as $\epsilon$ is varied we
obtain Fig. 11. As shown by the circles (from direct simulations of Eq. (39))
this solution is not always stable. As $\epsilon$ is decreased the solution
loses stability to a uniformly travelling wave, and the period of this wave is
not plotted.

Figure 11: Period, $T$, of periodic solutions of Eq. (39) of the form shown in
Fig. 10. Circles are from direct simulations of Eq. (39). Parameters:
$\omega_{0}=1$, $\Delta=0.1$.

##  IV Discussion

We considered time-periodic solutions of the Eqs. (4)–(5), which exactly
describe the asymptotic dynamics of the network (1) in the limit of
$N\to\infty$. At every point in space, Eq. (4) is a Riccati equation and we
used this to derive a self-consistency equation that every periodic solution
of Eq. (4) must satisfy. The Poincaré map of the Riccati equation is a Möbius
map and we can determine this map at every point in space using just three
numerical solutions of the Riccati equation. Knowing the Möbius map enables us
to numerically solve the self-consistency equation in a computationally
efficient manner. We showed the results of numerically continuing several
types of periodic solutions as a parameter was varied.

We derived equations governing the stability of such periodic solutions, but
solving these equations is numerically challenging. We also derived the
expression for the mean firing rate of neurons in a network in terms of the
quantities already calculated in the self-consistency equation. We finished in
Sec. III by demonstrating the application of our approach to several other
models involving delays, two populations of neurons, and a network of Winfree
oscillators.

Our approach relies critically on the mathematical form of the continuum-level
equations (they can be written as a Riccati equation) which are derived using
the Ott/Antonsen ansatz, valid only for phase oscillators whose dynamics and
coupling involve sinusoidal functions of phases or phase differences. Other
systems for which our approach should be applicable include two-dimensional
networks which support moving or “breathing” solutions [5]; however the
coupling function would have to be of the form such that the integral
equivalent to (5) could be written exactly using a small set of spatial basis
functions. Another application would be to any system which is periodically
forced in time and responds in a periodic way [54, 50, 53].

## Appendix

Let us consider a complex Riccati equation

$\frac{\displaystyle dz}{\displaystyle dt}=c_{0}(t)+c_{1}(t)z+c_{2}(t)z^{2}$
(42)

with $2\pi$-periodic complex-valued coefficients $c_{0}(t)$, $c_{1}(t)$ and
$c_{2}(t)$. In this section we prove a statement which is a modified version
of Proposition 3 from [43].

###### Proposition 2

Suppose that there is $c_{*}>0$ such that

$\mathrm{Re}\>(\overline{c_{0}(t)}z+c_{1}(t)+c_{2}(t)z)\leq-
c_{*}\mathrm{Re}\>(z+1)$ (43)

for all $z\in\overline{\mathbb{D}}$ and $0\leq t\leq 2\pi$, and that $z=-1$ is
not a fixed point of Eq. (42). Then, the Poincaré map of Eq. (42) is described
by a hyperbolic or loxodromic Möbius transformation. Moreover, the stable
fixed point of this map lies in the open unit disc $\mathbb{D}$, while the
unstable fixed point lies in the complementary domain
$\hat{\mathbb{C}}\backslash\overline{\mathbb{D}}$, where
$\hat{\mathbb{C}}=\mathbb{C}\cup\\{\infty\\}$ is the extended complex plane.
For Eq. (42) this means that it has exactly one stable $2\pi$-periodic
solution and this solution satisfies $|z(t)|<1$ for all $0\leq t\leq 2\pi$.

Proof: The fact that the Poincaré map of Eq. (42) is described by a Möbius
transformation was shown elsewhere [10, 59]. So we only need to show that such
a Möbius transformation $\mathcal{M}(z)$ maps all points of the unit circle
$|z|=1$ into the open unit disc $\mathbb{D}$ (but not on its boundary). Then
we can repeat the arguments of the proof of Proposition 2 from [43].

Suppose the opposite. Then Eq. (42) has a solution $z_{*}(t)$ such that
$|z_{*}(0)|=|z_{*}(2\pi)|=1$. Due to the inequality (43) and the Proposition 1
from [43] this solution satisfies $|z_{*}(t)|\leq 1$ for all $t\in[0,2\pi]$.
Moreover, since $z=-1$ is not a fixed point of Eq. (42), we can always choose
$t_{*}\in(0,2\pi)$ such that $z_{*}(t)\neq-1$ for $t\in[t_{*},2\pi)$. Then the
Mean Value Theorem implies

$0\leq|z_{*}(2\pi)|^{2}-|z_{*}(t_{*})|^{2}=2\pi\frac{\displaystyle
d|z_{*}|^{2}}{\displaystyle dt}(t_{**})$ (44)

for some $t_{**}\in(t_{*},2\pi)$. On the other hand, due to our assumptions,
we have

$\displaystyle\frac{\displaystyle d|z_{*}|^{2}}{\displaystyle dt}$
$\displaystyle=$ $\displaystyle
2\mathrm{Re}\>\left(\overline{c_{0}(t)}z_{*}(t)+c_{1}(t)+c_{2}(t)z_{*}(t)\right)$
$\displaystyle\leq$ $\displaystyle-2c_{*}\mathrm{Re}\>(z_{*}(t)+1)$

and therefore

$\frac{\displaystyle d|z_{*}|^{2}}{\displaystyle
dt}(t_{**})\leq-2c_{*}\mathrm{Re}\>(z_{*}(t_{**})+1)<0.$

This is a contradiction to (44), which completes the proof.

###### Remark 1

If the conditions of Proposition 2 are fulfilled, then the stable solution of
Eq. (42) is also asymptotically stable. This result follows from the
properties of stable fixed points of hyperbolic and loxodromic Möbius
transformations.

###### Remark 2

For every fixed $x$ the equation (7) from the main text is equivalent to Eq.
(42) with

$c_{0}(t)=c_{2}(t)=i\left(W(x,t)+\zeta-\frac{\displaystyle 1}{\displaystyle
2\omega}\right)$

and

$c_{1}(t)=2i\left(W(x,t)+\zeta+\frac{\displaystyle 1}{\displaystyle
2\omega}\right),$

where $W(x,t)$ is a real-valued function, $\omega$ is a positive constant, and
$\zeta=(\eta_{0}+i\gamma)/(2\omega)$ with $\eta_{0}\in\mathbb{R}$ and
$\gamma>0$. In this case, we have

$\mathrm{Re}\>(\overline{c_{0}(t)}z+c_{1}(t)+c_{2}(t)z)=-\frac{\displaystyle\gamma}{\displaystyle\omega}\mathrm{Re}\>(z+1)$

and therefore inequality (43) is satisfied. On the other hand, we have

$c_{0}(t)+c_{1}(t)(-1)+c_{2}(t)(-1)^{2}=-2i/\omega\neq 0$

and therefore $z=-1$ is not a fixed point of the corresponding Riccati
equation. Thus, all of the conclusions of Proposition 2 hold true for Eq. (7).

###### Remark 3

If the conditions of Proposition 2 are fulfilled, then the stable solution of
Eq. (42) can be computed in the following way.

(i) One solves Eq. (42) on the interval $t\in(0,2\pi]$ with three different
initial conditions $z(0)=z_{k}\in\mathbb{D}$, $k=1,2,3$, and obtains three
solutions $Z_{k}(t)$. Since each $z_{k}$ lies in the open unit disc
$\mathbb{D}$ this automatically implies $|Z_{k}(t)|<1$ for all $t\in(0,2\pi]$.

(ii) One denotes $w_{k}=Z_{k}(2\pi)$. Then, due to the properties of Poincaré
map one has $w_{k}=\mathcal{M}(z_{k})$, $k=1,2,3$, where $\mathcal{M}(z)$ is a
Möbius transformation representing this map. The above three relations can be
used to reconstruct the map $\mathcal{M}(z)$, namely

$\mathcal{M}(z)=\frac{az+b}{cz+d}$

where

$\displaystyle
a=\det\left(\begin{array}[]{ccc}z_{1}w_{1}&w_{1}&1\\\\[5.69054pt]
z_{2}w_{2}&w_{2}&1\\\\[5.69054pt] z_{3}w_{3}&w_{3}&1\end{array}\right),\quad
b=\det\left(\begin{array}[]{ccc}z_{1}w_{1}&z_{1}&w_{1}\\\\[5.69054pt]
z_{2}w_{2}&z_{2}&w_{2}\\\\[5.69054pt]
z_{3}w_{3}&z_{3}&w_{3}\end{array}\right),$ $\displaystyle
c=\det\left(\begin{array}[]{ccc}z_{1}&w_{1}&1\\\\[5.69054pt]
z_{2}&w_{2}&1\\\\[5.69054pt]
z_{3}&w_{3}&1\end{array}\right),\quad\phantom{w_{1}}d=\det\left(\begin{array}[]{ccc}z_{1}w_{1}&z_{1}&1\\\\[5.69054pt]
z_{2}w_{2}&z_{2}&1\\\\[5.69054pt] z_{3}w_{3}&z_{3}&1\end{array}\right).$

(iii) Once the map $\mathcal{M}(z)$ is known, one can find its fixed points by
solving the quadratic equation

$cz^{2}+dz-az-b=0.$

This yields two roots

$z_{\pm}=\frac{\displaystyle a-d\pm\sqrt{(a-d)^{2}+4bc}}{\displaystyle 2c}.$

(iv) Choosing from the roots $z_{+}$ and $z_{-}$ the one that lies in the unit
disc $\mathbb{D}$, one obtains the initial condition that determines the
periodic solution of interest. The latter can be computed by solving Eq. (42)
with this initial condition.

(v) Sometime it may happen that the Poincaré map $\mathcal{M}(z)$ is strongly
contracting so that

$|w_{1}-w_{2}|+|w_{3}-w_{2}|<10^{-8},$

where the value $10^{-8}$ is chosen through experience. In this case, the
calculations in steps (ii) and (iii) become inaccurate. Then the initial
condition of the periodic solution of interest is approximately given by the
average $(w_{1}+w_{2}+w_{3})/3$.

The above steps (i)–(v) can be understood as a constructive definition of the
solution operator of Eq. (42), which for every admissible choise of
$2\pi$-periodic coefficients $c_{1}(t)$, $c_{2}(t)$ and $c_{3}(t)$ yields the
corresponding stable $2\pi$-periodic solution of Eq. (42). More detailed
justification of this definition can be found in [43, Section 4]. The
algorithm in Sec. II.1 consists of applying the procedure above at every point
on the spatial grid (in parallel).

## Acknowledgements

The work of O.E.O. was supported by the Deutsche Forschungsgemeinschaft under
Grant No. OM 99/2-2.

## Declarations

Competing interests: The authors have no competing interests to declare that
are relevant to the content of this article.

Code availability: Code for calculating any of the results presented here is
available on reasonable request from the authors.

## Author contributions

#### Conceptualization:

C.R. Laing, O.E. Omel’chenko.

#### Formal analysis:

C.R. Laing, O.E. Omel’chenko.

#### Funding acquisition:

O.E. Omel’chenko.

#### Investigation:

C.R. Laing, O.E. Omel’chenko.

#### Methodology:

C.R. Laing, O.E. Omel’chenko.

#### Software:

C.R. Laing.

#### Validation:

C.R. Laing, O.E. Omel’chenko.

#### Visualization:

C.R. Laing.

#### Writing – original draft:

C.R. Laing, O.E. Omel’chenko.

#### Writing – review & editing:

C.R. Laing, O.E. Omel’chenko.

## References

  * [1] Amari, S.i.: Dynamics of pattern formation in lateral-inhibition type neural fields.  Biol. Cybern. 27(2), 77–87 (1977)
  * [2] Ariaratnam, J.T., Strogatz, S.H.: Phase diagram for the winfree model of coupled nonlinear oscillators.  Phys. Rev. Lett. 86(19), 4278 (2001)
  * [3] Atay, F.M., Hutt, A.: Neural fields with distributed transmission speeds and long-range feedback delays.  SIAM J. Appl. Dyn. Syst. 5(4), 670–698 (2006)
  * [4] Avitabile, D., Desroches, M., Ermentrout, G.B.: Cross-scale excitability in networks of quadratic integrate-and-fire neurons.  PLoS Comput. Biol. 18(10), e1010569 (2022)
  * [5] Bataille-Gonzalez, M., Clerc, M.G., Omel’chenko, O.E.: Moving spiral wave chimeras.  Phys. Rev. E 104, L022203 (2021)
  * [6] Blomquist, P., Wyller, J., Einevoll, G.T.: Localized activity patterns in two-population neuronal networks.  Phys. D: Nonlinear Phenom. 206(3), 180–212 (2005)
  * [7] Börgers, C., Kopell, N.: Effects of noisy drive on rhythms in networks of excitatory and inhibitory neurons.  Neural Comput. 17(3), 557–608 (2005)
  * [8] Byrne, A., Avitabile, D., Coombes, S.: Next-generation neural field model: The evolution of synchrony within patterns and waves.  Phys. Rev. E 99, 012313 (2019)
  * [9] Byrne, Á., Ross, J., Nicks, R., Coombes, S.: Mean-field models for EEG/MEG: from oscillations to waves.  Brain Topogr. 35(1), 36–53 (2022)
  * [10] Campos, J.: Möbius transformations and periodic solutions of complex Riccati equations.  Bull. London Math. Soc. 29, 205–215 (1997)
  * [11] Coombes, S., Byrne, Á.: Next generation neural mass models.  In: Nonlinear Dynamics in Computational Neuroscience, pp. 1–16. Springer (2019)
  * [12] Coombes, S., Laing, C.: Delays in activity-based neural networks.  Philos. Trans. Royal Soc. A: Math. Phys. Eng. Sci. 367(1891), 1117–1129 (2009)
  * [13] Daleckiĭ, J.L., Kreĭn, M.G.: Stability of solutions of differential equations in Banach space.  43\. American Mathematical Soc. (2002)
  * [14] Devalle, F., Montbrió, E., Pazó, D.: Dynamics of a large system of spiking neurons with synaptic delay.  Phys. Rev. E 98(4), 042214 (2018)
  * [15] Devalle, F., Roxin, A., Montbrió, E.: Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks.  PLoS Comput. Biol. 13(12), e1005881 (2017)
  * [16] Ermentrout, G.B., Kopell, N.: Parabolic bursting in an excitable system coupled with a slow oscillation.  SIAM J. Appl. Math. 46(2), 233–253 (1986)
  * [17] Esnaola-Acebes, J.M., Roxin, A., Avitabile, D., Montbrió, E.: Synchrony-induced modes of oscillation of a neural field model.  Phys. Rev. E 96, 052407 (2017)
  * [18] Funahashi, S., Bruce, C.J., Goldman-Rakic, P.S.: Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex.  J. Neurophysiol. 61(2), 331–349 (1989)
  * [19] Gallego, R., Montbrió, E., Pazó, D.: Synchronization scenarios in the Winfree model of coupled oscillators.  Phys. Rev. E 96(4), 042208 (2017)
  * [20] Jiruska, P., de Curtis, M., Jefferys, J.G.R., Schevon, C.A., Schiff, S.J., Schindler, K.: Synchronization and desynchronization in epilepsy: controversies and hypotheses.  J. Physiol. 591(4), 787–797 (2013)
  * [21] Keane, A., Krauskopf, B., Dijkstra, H.A.: The effect of state dependence in a delay differential equation model for the el niño southern oscillation.  Philos. Trans. Royal Soc. A: Math. Phys. Eng. Sci. 377(2153), 20180121 (2019)
  * [22] Kilpatrick, Z.P., Ermentrout, B.: Wandering bumps in stochastic neural fields.  SIAM J. Appl. Dyn. Syst. 12(1), 61–94 (2013)
  * [23] Kuramoto, Y.: Chemical Oscillations, Waves, and Turbulence.  Springer, Berlin (1984)
  * [24] Laing, C., Chow, C.: Stationary bumps in networks of spiking neurons.  Neural Comput. 13(7), 1473–1494 (2001)
  * [25] Laing, C.R.: Derivation of a neural field model from a network of theta neurons.  Phys. Rev. E 90(1), 010901 (2014)
  * [26] Laing, C.R.: Numerical bifurcation theory for high-dimensional neural models.  J. Math. Neurosci. 4(1), 13 (2014)
  * [27] Laing, C.R.: Exact neural fields incorporating gap junctions.  SIAM J. Appl. Dyn. Syst. 14(4), 1899–1929 (2015)
  * [28] Laing, C.R.: Phase oscillator network models of brain dynamics.  In: A. Moustafa (ed.) Computational Models of Brain and Behavior, chap. 37, pp. 505–517. Wiley-Blackwell, Hoboken, NJ (2017)
  * [29] Laing, C.R., Bläsche, C., Means, S.: Dynamics of structured networks of Winfree oscillators.  Front. Syst. Neurosci. 15, 631377 (2021)
  * [30] Laing, C.R., Longtin, A.: Dynamics of deterministic and stochastic paired excitatory-inhibitory delayed feedback.  Neural Comput. 15(12), 2779–2822 (2003)
  * [31] Laing, C.R., Omel’chenko, O.: Moving bumps in theta neuron networks.  Chaos 30(4), 043117 (2020)
  * [32] Laing, C.R., Troy, W.: PDE methods for nonlocal models.  SIAM J. Appl. Dyn. Syst. 2(3), 487–516 (2003)
  * [33] Laing, C.R., Troy, W.C., Gutkin, B., Ermentrout, G.B.: Multiple bumps in a neuronal model of working memory.  SIAM J. Appl. Math. 63(1), 62–97 (2002)
  * [34] Latham, P., Richmond, B., Nelson, P., Nirenberg, S.: Intrinsic dynamics in neuronal networks. I. Theory.  J. Neurophysiol. 83(2), 808–827 (2000)
  * [35] Lee, W.S., Ott, E., Antonsen, T.M.: Large coupled oscillator systems with heterogeneous interaction delays.  Phys. Rev. Lett. 103, 044101 (2009)
  * [36] Lindén, H., Petersen, P.C., Vestergaard, M., Berg, R.W.: Movement is governed by rotational neural dynamics in spinal motor networks.  Nature 610, 526–531 (2022)
  * [37] Montbrió, E., Pazó, D., Roxin, A.: Macroscopic description for networks of spiking neurons.  Phys. Rev. X 5, 021028 (2015)
  * [38] Netoff, T.I., Schiff, S.J.: Decreased neuronal synchronization during experimental seizures.  J. Neurosci. 22(16), 7297–7307 (2002)
  * [39] Omel’chenko, O.: Mathematical framework for breathing chimera states.  J. Nonlinear Sci. 32(2), 1–34 (2022)
  * [40] Omel’chenko, O., Laing, C.R.: Collective states in a ring network of theta neurons.  Proc. Royal Soc. A 478(2259), 20210817 (2022)
  * [41] Omel’chenko, O., Wolfrum, M., Laing, C.R.: Partially coherent twisted states in arrays of coupled phase oscillators.  Chaos 24, 023102 (2014)
  * [42] Omel’chenko, O.E.: Nonstationary coherence-incoherence patterns in nonlocally coupled heterogeneous phase oscillators.  Chaos 30(4), 043103 (2020)
  * [43] Omel’chenko, O.E.: Periodic orbits in the Ott-Antonsen manifold.  Nonlinearity 36, 845–861 (2023)
  * [44] Ott, E., Antonsen, T.M.: Low dimensional behavior of large systems of globally coupled oscillators.  Chaos 18(3), 037113 (2008)
  * [45] Ott, E., Antonsen, T.M.: Long time evolution of phase oscillator systems.  Chaos 19(2), 023117 (2009)
  * [46] Pazó, D., Montbrió, E.: Low-dimensional dynamics of populations of pulse-coupled oscillators.  Phys. Rev. X 4, 011009 (2014)
  * [47] Pietras, B.: Pulse shape and voltage-dependent synchronization in spiking neuron networks.  arXiv preprint arXiv:2304.09813 (2023)
  * [48] Pinto, D.J., Ermentrout, G.B.: Spatially structured activity in synaptically coupled neuronal networks: II. Lateral inhibition and standing pulses.  SIAM J. Appl. Math. 62(1), 226–243 (2001)
  * [49] Ratas, I., Pyragas, K.: Macroscopic self-oscillations and aging transition in a network of synaptically coupled quadratic integrate-and-fire neurons.  Phys. Rev. E 94, 032215 (2016)
  * [50] Reyner-Parra, D., Huguet, G.: Phase-locking patterns underlying effective communication in exact firing rate models of neural networks.  PLOS Comput. Biol. 18(5), e1009342 (2022)
  * [51] Roxin, A., Brunel, N., Hansel, D.: Role of delays in shaping spatiotemporal dynamics of neuronal activity in large networks.  Phys. Rev. Lett. 94(23), 238103 (2005)
  * [52] Schmidt, H., Avitabile, D.: Bumps and oscillons in networks of spiking neurons.  Chaos 30, 033133 (2020)
  * [53] Schmidt, H., Avitabile, D., Montbrió, E., Roxin, A.: Network mechanisms underlying the role of oscillations in cognitive tasks.  PLoS Comput. Biol. 14(9), e1006430 (2018)
  * [54] Segneri, M., Bi, H., Olmi, S., Torcini, A.: Theta-nested gamma oscillations in next generation neural mass models.  Front. Comput. Neurosci. 14, 47 (2020)
  * [55] Shima, S., Kuramoto, Y.: Rotating spiral waves with phase-randomized core in nonlocally coupled oscillators.  Phys. Rev. E 69(3), 036213 (2004)
  * [56] Strogatz, S.: From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators.  Phys. D 143(1-4), 1–20 (2000)
  * [57] Uhlhaas, P.J., Singer, W.: Abnormal neural oscillations and synchrony in schizophrenia.  Nature Rev. Neurosci. 11(2), 100–113 (2010)
  * [58] di Volo, M., Torcini, A.: Transition from asynchronous to oscillatory dynamics in balanced spiking networks with instantaneous synapses.  Phys. Rev. Lett. 121, 128301 (2018)
  * [59] Wilczyński, P.: Planar nonautonomous polynomial equations: the riccati equation.  J. Differ. Equ. 244, 1304–1328 (2008)
  * [60] Wimmer, K., Nykamp, D.Q., Constantinidis, C., Compte, A.: Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory.  Nature Neurosci. 17(3), 431–439 (2014)
  * [61] Zhang, K.: Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory.  J. Neurosci. 16(6), 2112–2126 (1996)