# Activity patterns in ring networks of quadratic integrate-and-fire neurons
with synaptic and gap junction coupling

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

###### Abstract

We consider a ring network of quadratic integrate-and-fire neurons with
nonlocal synaptic and gap junction coupling. The corresponding neural field
model supports solutions such as standing and travelling waves, and also
lurching waves. We show that many of these solutions satisfy self-consistency
equations which can be used to follow them as parameters are varied. We
perform numerical bifurcation analysis of the neural field model,
concentrating on the effects of varying gap junction coupling strength. Our
methods are generally applicable to a wide variety of networks of quadratic
integrate-and-fire neurons.

neural field model, Riccati equation, quadratic integrate-and-fire neuron,
Ott/Antonsen, self-consistency

##  I Introduction

The collective behaviour of spatially extended networks of neurons is a topic
of ongoing interest [16, 13, 12, 3, 7, 14]. While it is possible to simulate
large networks of model neurons [44, 22], continuum level descriptions (neural
field models) often provide more potential for mathematical analysis. Early
neural field models were phenomenological [1, 30, 43, 11] but more recently
models derived rigorously from infinite networks of spiking neurons have
become available [27, 36, 25, 29, 17, 40, 6].

One type of solution of a neural field model is a “bump” — a spatially
localised group of active neurons. Such solutions are thought to be relevant
to working memory [9] and the head direction system [45, 24]. Also important
are travelling waves [38, 5, 23, 15]. These are thought to be relevant for
modelling epilepsy and migraines, see for example [4]. We are often interested
in when and how such solutions lose stability or are destroyed in
bifurcations, and what patterns are stable beyond such bifurcations. For
example, a bump may start to “breathe” [18] or a travelling wave may no longer
have a fixed profile.

In this paper we study a ring network of quadratic integrate-and-fire (QIF)
neurons coupled both synaptically and via gap junctions. We analyse the
continuum network, whose dynamics are described by a neural field model. In
similar previous work we considered networks of theta neurons with just
synaptic coupling, concentrating on only time-periodic solutions [29]. (The
QIF neuron with infinite threshold and reset is equivalent to a theta neuron
[33].) In earlier work we also studied the stationary states of a ring network
of synaptically coupled theta neurons [36]. On the other hand, models similar
(but not identical) to ours have been the focus of other reasearchers. Byrne
et al. studied networks of QIF neurons with both synaptic and gap junction
coupling, but their synaptic coupling incorporated propagation delays and the
gap junction coupling was purely local [7]. Their neural field model was only
valid in the long wavelength limit. Schmidt and Avitabile considered a ring
network of QIF neurons with nonlocal synaptic coupling and analysed both
steady states and time dependent solutions which arose as the result of
periodic forcing [40]. Esnaola-Acebes et al. studied a similar model,
focussing on the decaying oscillatory modes that such models show [17]. Byrne
et al. also considered a ring network of QIF neurons with synaptic coupling,
performing a largely numerical study of some of the possible types of
solutions [6].

Our model is probably most similar to that of [7], but one of our
contributions is to show that some solutions of interest can be analysed using
a self-consistency approach that has previously been used only for networks of
theta neurons or Kuramoto and Winfree oscillators [29, 35, 2]. The structure
of the paper is as follows. In Sec. II we present the model and describe the
types of patterns typically seen when numerically solving it. In Sec. III we
choose two values of the synaptic coupling strength and vary the strength of
the gap junction coupling, observing transitions between various types of
solutions. In Sec. IV we present analytical methods which can be used to
describe most of the solutions observed in Sec. III. Section V shows the
results of implementing the methods described in Sec. IV and includes a
discussion of the relative efficiencies of several alternative methods. We
conclude in Sec. VI, and the Appendix contains a number of useful results
about the solutions of the complex Riccati equation.

##  II Model and types of solutions

Let us consider a spatially extended network of $N$ quadratic integrate-and-
fire (QIF) neurons with both gap-junction and pulsatile synaptic coupling:

$\frac{\displaystyle dV_{j}}{\displaystyle
dt}=\eta_{j}+V_{j}^{2}+\kappa_{\mathrm{v}}\ell\sum\limits_{k=1}^{N}W_{\mathrm{v}}\left(|j-k|\ell\right)(V_{k}-V_{j})+\kappa_{\mathrm{s}}\ell\sum\limits_{k=1}^{N}\sum\limits_{m\in\mathbb{Z}}W_{\mathrm{s}}\left(|j-k|\ell\right)\delta\left(t-T_{k}^{m}\right).$
(1)

Here, the spiking events of the $j$th neuron $T_{j}^{m}$ are determined by the
reset condition: if $V_{j}\to+\infty$ for $t\nearrow T_{j}^{m}$, then
$V_{j}\to-\infty$ for $t\searrow T_{j}^{m}$.

We assume that the network is organized as a one-dimensional array with
periodic boundary conditions and distance-dependent connectivity defined by
coupling functions $W_{\mathrm{v}}$ and $W_{\mathrm{s}}$. The neurons differ
from each other only in the excitability parameters $\eta_{j}$, which are
chosen randomly and independently from a Lorentzian distribution

$g(\eta)=\frac{\displaystyle 1}{\displaystyle\pi}\frac{\displaystyle
1}{\displaystyle(\eta-\eta_{0})^{2}+\gamma^{2}}\quad\mbox{with}\quad\eta_{0}\in\mathbb{R}\quad\mbox{and}\quad\gamma>0.$

The coupling functions $W_{\mathrm{v}}$ and $W_{\mathrm{s}}$ are assumed to be
even, continuous and periodic

$W_{\mathrm{v}}\left((k+N)\ell\right)=W_{\mathrm{v}}\left(k\ell\right),\qquad
W_{\mathrm{s}}\left((k+N)\ell\right)=W_{\mathrm{s}}\left(k\ell\right),$

and the intensity of interactions between neurons is controlled by two scalar
coupling strengths $\kappa_{\mathrm{v}}$ and $\kappa_{\mathrm{s}}$. Moreover,
for the sake of simplicity of the analytical consideration below, it is
convenient to assume $\ell=2\pi/N$. In this case, functions
$W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$ are $2\pi$-periodic and can be
represented in the form of Fourier series

$W_{\mathrm{v}}(x)=\sum\limits_{m=-\infty}^{\infty}\hat{W}_{\mathrm{v},m}e^{imx}\quad\mbox{and}\quad
W_{\mathrm{s}}(x)=\sum\limits_{m=-\infty}^{\infty}\hat{W}_{\mathrm{s},m}e^{imx}$

with coefficients

$\hat{W}_{\mathrm{v},m}=\frac{\displaystyle 1}{\displaystyle
2\pi}\int_{-\pi}^{\pi}W_{\mathrm{v}}(x)e^{-imx}dx\quad\mbox{and}\quad\hat{W}_{\mathrm{s},m}=\frac{\displaystyle
1}{\displaystyle 2\pi}\int_{-\pi}^{\pi}W_{\mathrm{s}}(x)e^{-imx}dx.$ (2)

Note that for even functions $W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$, the
coefficients $\hat{W}_{\mathrm{v},m}$ and $\hat{W}_{\mathrm{s},m}$ are real
and therefore

$\hat{W}_{\mathrm{v},m}=\frac{\displaystyle 1}{\displaystyle
2\pi}\int_{-\pi}^{\pi}W_{\mathrm{v}}(x)\cos(mx)dx\quad\mbox{and}\quad\hat{W}_{\mathrm{s},m}=\frac{\displaystyle
1}{\displaystyle 2\pi}\int_{-\pi}^{\pi}W_{\mathrm{s}}(x)\cos(mx)dx.$

It is well-known [27, 7, 37, 32] that in the continuum limit $N\to\infty$, the
long-term dynamics of system (1) can be described by a neural field equation

$\frac{\displaystyle\partial u}{\displaystyle\partial
t}=\gamma-\kappa_{\mathrm{v}}u+i\left[\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(u)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(u)-u^{2}\right]$
(3)

for a complex-valued function $u=u(x,t)$, where $\mathcal{K}_{\mathrm{v}}$ and
$\mathcal{K}_{\mathrm{s}}$ denote two integral operators

$(\mathcal{K}_{\mathrm{v}}\varphi)(x)=\int_{0}^{2\pi}W_{\mathrm{v}}(x-y)\varphi(y)dy\quad\mbox{and}\quad(\mathcal{K}_{\mathrm{s}}\varphi)(x)=\int_{0}^{2\pi}W_{\mathrm{s}}(x-y)\varphi(y)dy.$

Note that the real and imaginary parts of the solution $u(x,t)$ can be related
to the local firing rate

$R(x,t)=\frac{\displaystyle 1}{\displaystyle\pi}\mathrm{Re}\>u(x,t)$ (4)

and the local mean field potential

$V(x,t)=\mathrm{Im}\>u(x,t).$ (5)

Due to this interpretation, only those solutions that satisfy
$\mathrm{Re}\>(u)\geq 0$ are physically meaningful. Moreover, it can be
explicitly shown that the property $\mathrm{Re}\>(u)\geq 0$ is preserved by
the dynamics of Eq. (3). Therefore, we focus on this invariant set below.

Note also that the mapping

$z(x,t)=\frac{\displaystyle 1-\overline{u}}{\displaystyle 1+\overline{u}},$
(6)

where overline denotes the complex conjugate, transforms the function $u$ into
an analog of the Kuramoto local order parameter, $z$; the dynamics can be
written equally-well in terms of $z$ [27, 7].

In the examples below, we choose $W_{\mathrm{v}}(x)$ in the form of a Gaussian
function

$W_{\mathrm{v}}(x)=\frac{\displaystyle
1}{\displaystyle\sqrt{2\pi}\sigma}e^{-x^{2}/(2\sigma^{2})}\quad\mbox{with}\quad\sigma=0.1,$
(7)

and $W_{\mathrm{s}}(x)$ in the form of a Mexican-hat function which is
positive for small $|x|$ and negative for larger $|x|$:

$W_{\mathrm{s}}(x)=\frac{\displaystyle
1}{\displaystyle\sqrt{2\pi}\sigma_{1}}e^{-x^{2}/(2\sigma_{1}^{2})}-\frac{\displaystyle
1}{\displaystyle\sqrt{2\pi}\sigma_{2}}e^{-x^{2}/(2\sigma_{2}^{2})}\quad\mbox{with}\quad\sigma_{1}=0.5,\>\sigma_{2}=1.$
(8)

More precisely, the above functions determine the coupling only for
$|x|\leq\pi$, while we use their $2\pi$-periodic extensions for other values
of $x$. Note that the width of the gap junction coupling kernel,
$W_{\mathrm{v}}$, is much smaller than that of the synaptic coupling kernel,
$W_{\mathrm{s}}$, since gap junctional coupling is typically more localised
than synaptic coupling [27].

In numerical simulations for the spatially discretized version of Eq. (3), we
usually encountered the following five types of stable states:

  1. 1.

Uniform states, i.e. constant solutions of Eq. (3) that do not depend on
either $x$ or $t$.

  2. 2.

Stationary states, i.e. time-independent solutions of the form $u=a(x)$ for
some $2\pi$-periodic function $a$.

  3. 3.

Traveling waves, i.e. solutions of the form $u=a(x-st)$, which are stationary
in a frame moving with a constant speed $s\neq 0$.

  4. 4.

Standing waves, i.e. solutions $u=a(x,t)$ satisfying the periodicity condition
$a(x,t+T)=a(x,t)$ with some $T>0$.

  5. 5.

Lurching waves, i.e. solutions $u=a(x,t)$ satisfying the shifted periodicity
condition $a(x+\chi,t+T)=a(x,t)$ with some $\chi\neq 0$ and $T>0$. Such waves
have been found previously in neural models [41, 10, 19, 20, 21, 46].

In this paper we analyse all of these types of solutions using a mixture of
numerical and analytical methods. In Sec. III we show examples of these types
of solution and characterise some aspects of them as the parameter
$\kappa_{\mathrm{v}}$ (the strength of gap junction coupling) is varied. We
choose to vary this parameter as the influence of gap junction coupling in
neural field models has only recently been considered [27, 32, 37]. In Sec. IV
we undertake bifurcation analysis of the various types of solution, sometimes
describing them using a recently presented method that can be used to
efficiently characterise solutions using self-consistency arguments [35, 34,
29, 2].

##  III Phenomenology

To determine the typical attractors of the neural field equation (3), we
discretized it on a uniform spatial grid of $1024$ nodes and performed
numerical simulations for various parameters $\kappa_{\mathrm{s}}$ and
$\kappa_{\mathrm{v}}$. Other parameters were chosen as $\eta_{0}=1$ and
$\gamma=0.5$, by analogy with a similar model considered by Byrne et al. in
[7]. Our observations are summarized below, where we provide a comprehensive
overview of two representative cases $\kappa_{\mathrm{s}}=10$ and
$\kappa_{\mathrm{s}}=20$.

###  III.1 Case $\kappa_{\mathrm{s}}=10$

Choosing $\eta_{0}=1$, $\gamma=0.5$, $\kappa_{\mathrm{s}}=10$ and
$\kappa_{\mathrm{v}}=1$ and running the numerical simulations for Eq. (3) with
various initial conditions, we found two coexisting spatio-temporal patterns;
see Fig. 1. One of them had the appearance of a periodically oscillating
standing wave, while the other had the appearance of a travelling wave,
although the wave does not travel with a constant profile. Note that standing
and travelling waves have been observed in several other neural field models
[29, 6, 7].

Figure 1: (a) A standing wave solution of Eq. (3) and (b) a coexisting
travelling wave solution. For both states, the local firing rate $R$,
calculated by formula (4), is shown. Parameters: $\kappa_{\mathrm{s}}=10$,
$\kappa_{\mathrm{v}}=1$, $\eta_{0}=1$ and $\gamma=0.5$.

Next, by adiabatically changing the coupling strength $\kappa_{\mathrm{v}}$,
we investigated the stability intervals of these patterns. (When moving from
one value of $\kappa_{\mathrm{v}}$ to the next, we added a small random
perturbation to the final state to give the new initial state.) To this end,
for each value of $\kappa_{\mathrm{v}}$, we calculated a trajectory of length
$1000$ time units, discarding the previous transient of the same length. To
characterize the dynamics of the obtained trajectory, we used the spatially
averaged firing rate

$\langle R(t)\rangle=\frac{\displaystyle 1}{\displaystyle
2\pi}\int_{0}^{2\pi}R(x,t)dx=\frac{\displaystyle 1}{\displaystyle
2\pi^{2}}\int_{0}^{2\pi}\mathrm{Re}\>u(x,t)dx$

and its time variation

$\Delta R=\max\limits_{t}\langle R(t)\rangle-\min\limits_{t}\langle
R(t)\rangle.$

where $t$ runs over the length of the simulation. In the case $\Delta R>0$,
i.e. when $\langle R(t)\rangle$ is non-constant, we used the plot of $\langle
R(t)\rangle$ to calculate the time intervals $\Delta t$ between consecutive
local maxima of this function. The resulting dependences of $\Delta R$ and
$\Delta t$ on $\kappa_{\mathrm{v}}$ are shown in Fig. 2. Moreover, for every
time-periodic solution of Eq. (3), we also show its least period $T$, which
was calculated by averaging the time intervals between the consecutive local
maxima of $\mathrm{Re}\>u(0,t)$. (Note that for travelling waves we needed to
double the above value, since $\mathrm{Re}\>u(0,t)$ attains two maxima during
the period. Moreover, for travelling waves with $\kappa_{\mathrm{v}}>0.964$ we
also took into account that additional local maxima with $\mathrm{Re}\>\
u(0,t)<1.5$ occur between the primary local maxima of $\mathrm{Re}\>\
u(0,t)$.)

Figure 2: Parameter sweeps of standing wave (a–c) and travelling wave (d–f)
solutions of Eq. (3). Three rows show the time variation of spatially averaged
firing rate $\Delta R$, the time intervals $\Delta t$ between the consecutive
local maxima in the plot of $\langle R(t)\rangle$, and the least period $T$.
Other parameters: $\eta_{0}=1$, $\gamma=0.5$ and $\kappa_{\mathrm{s}}=10$.

The diagrams in the left column of Fig. 2 indicate the existence of two types
of standing waves: the waves with a single maximal value of $\langle
R(t)\rangle$ and the waves with two alternating maximal values of $\langle
R(t)\rangle$. Typical examples of these patterns are shown in Fig. 3.

Figure 3: Two types of standing wave solutions of Eq. (3) for
$\kappa_{\mathrm{v}}=0.96$ (a) and $\kappa_{\mathrm{v}}=0.99$ (b). Other
parameters: $\eta_{0}=1$, $\gamma=0.5$ and $\kappa_{\mathrm{s}}=10$.

The diagrams in the right column of Fig. 2 have a more complicated structure.
There, we can distinguish four different types of travelling waves, see Fig.
4. The left-most values of $\kappa_{\mathrm{v}}$ correspond to a wave that
moves rigidly (notice $\Delta R=0$) with a constant speed (Fig. 4(a)). For a
larger value of $\kappa_{\mathrm{v}}$, we found a wave propagating with a non-
constant profile, which slowly drifts (Fig. 4(b)). This is a lurching wave.
For a still larger value of $\kappa_{\mathrm{v}}$ the lurching wave undergoes
some form of quasiperiodic modulation (Fig. 4(c)). Finally, for the right-most
values of $\kappa_{\mathrm{v}}$, we observe a wave moving above a turbulent
background (Fig. 4(d)).

Figure 4: Four types of travelling waves in Eq. (3) for
$\kappa_{\mathrm{v}}=0.96$ (a), $\kappa_{\mathrm{v}}=1.02$ (b),
$\kappa_{\mathrm{v}}=1.04$ (c) and $\kappa_{\mathrm{v}}=1.07$ (d). Other
parameters: $\eta_{0}=1$, $\gamma=0.5$ and $\kappa_{\mathrm{s}}=10$.

###  III.2 Case $\kappa_{\mathrm{s}}=20$

Keeping $\eta_{0}=1$, $\gamma=0.5$ and setting $\kappa_{\mathrm{s}}=20$ we
obtain several new types of solution for different values of
$\kappa_{\mathrm{v}}$; see Fig. 5. We see “two-bump” stationary solutions,
Fig. 5(a), for which there are two disjoint regions where neurons are firing
at a significant rate [31, 30]. We also see standing waves (or “breathing”
two-bump solutions, similar to those seen in [27, 6, 40]), Fig. 5(b), more
complex standing wave solutions, Fig. 5(c), lurching waves, Fig. 5(d), and
waves that travel but are not periodic in a uniformly-travelling coordinate
frame, Fig. 5(e,f). Quasi-statically sweeping through $\kappa_{\mathrm{v}}$ we
obtain Fig. 6.

Figure 5: Six types of solutions of Eq. (3) for $\kappa_{\mathrm{v}}=0.8$ (a),
$\kappa_{\mathrm{v}}=1$ (b), $\kappa_{\mathrm{v}}=1.3$ (c),
$\kappa_{\mathrm{v}}=1.6$ (d), $\kappa_{\mathrm{v}}=1.65$ (e) and
$\kappa_{\mathrm{v}}=1.8$ (f). Other parameters: $\eta_{0}=1$, $\gamma=0.5$
and $\kappa_{\mathrm{s}}=20$.

Figure 6: A parameter sweep for patterned solutions shown in Fig. 5. The
panels (a) and (b) show the time variation of spatially averaged firing rate
$\Delta R$ and the time intervals $\Delta t$ between the consecutive local
maxima in the plot of $\langle R(t)\rangle$, respectively. Other parameters:
$\eta_{0}=1$, $\gamma=0.5$ and $\kappa_{\mathrm{s}}=20$.

In the rest of the paper we analyse the types of solution shown, and provide a
more rigorous explanation for the results observed above.

##  IV Analytical methods

The kaleidoscopic set of spatiotemporal patterns found in Eq. (3) for various
system parameters can be logically explained by a bifurcation analysis of this
equation, focusing on its equilibrium and periodic solutions. We describe two
ways to do this. Our first approach relies on a discretized version of the
integro-differential equation (3) and uses standard software for the
bifurcation analysis of finite-dimensional dynamical systems. In the second
approach, we derive a set of integral self-consistency equations for different
types of solutions of Eq. (3) and carry out their finite-dimensional reduction
by Galerkin method.

###  IV.1 Uniform states

We start with the simplest form of solution, a uniform state. Recall that
these are constant solutions of Eq. (3), which do not depend on either $x$ or
$t$. On the other hand, these states will not be steady states of system (1),
as that describes a finite network of heterogeneous spiking neurons. However,
in system (1), such a state would have no macroscopic spatial structure nor
significant temporal variations.

####  IV.1.1 Existence

Each uniform state of Eq. (3) corresponds to a solution $u\in\mathbb{C}$ of
the equation

$\gamma-\kappa_{\mathrm{v}}u+i\left[\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(u)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(u)-u^{2}\right]=0.$
(9)

In the next proposition, we show that all physically meaningful solutions of
Eq. (9) lie on a certain one-parameter manifold.

###### Proposition 1

All constant solutions of Eq. (9) satisfying $\mathrm{Re}\>(u)\geq 0$ lie on
the manifold defined by formulas

$\displaystyle u$ $\displaystyle=$ $\displaystyle\frac{\displaystyle
1}{\displaystyle 2}f_{-}(F,\kappa_{\mathrm{v}},\gamma)+\frac{\displaystyle
i}{\displaystyle
2}\left(\kappa_{\mathrm{v}}-f_{+}(F,\kappa_{\mathrm{v}},\gamma)\right),$ (10)
$\displaystyle\eta_{0}$ $\displaystyle=$ $\displaystyle
F-\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},0}\left(\kappa_{\mathrm{v}}-f_{+}(F,\kappa_{\mathrm{v}},\gamma)\right)-\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},0}f_{-}(F,\kappa_{\mathrm{v}},\gamma),$
(11)

where

$f_{\pm}(F,\kappa_{\mathrm{v}},\gamma)=\sqrt{\frac{\displaystyle\sqrt{(\kappa_{\mathrm{v}}^{2}-4F)^{2}+16\gamma^{2}}\pm(\kappa_{\mathrm{v}}^{2}-4F)}{\displaystyle
2}},$

the coefficients $\hat{W}_{\mathrm{v},0}$ and $\hat{W}_{\mathrm{s},0}$ are
defined by the formula (2) with $m=0$, and $F\in\mathbb{R}$ is a free
parameter.

Proof: Let us rewrite Eq. (9) in the form

$\gamma+iF-\kappa_{\mathrm{v}}u-iu^{2}=0,$ (12)

where

$F=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(u)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(u).$
(13)

For a given $F$, Eq. (12) has one and only one solution satisfying
$\mathrm{Re}\>(u)\geq 0$, see Proposition 50. It is determined by the formula
(10). Inserting this result into Eq. (13) we obtain

$\displaystyle F$ $\displaystyle=$
$\displaystyle\eta_{0}+\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},0}\left(\kappa_{\mathrm{v}}-\sqrt{\frac{\displaystyle\sqrt{(\kappa_{\mathrm{v}}^{2}-4F)^{2}+16\gamma^{2}}+\kappa_{\mathrm{v}}^{2}-4F}{\displaystyle
2}}\right)$ $\displaystyle+$
$\displaystyle\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},0}\sqrt{\frac{\displaystyle\sqrt{(\kappa_{\mathrm{v}}^{2}-4F)^{2}+16\gamma^{2}}-\kappa_{\mathrm{v}}^{2}+4F}{\displaystyle
2}},$

what is an equivalent form of (11).

Note that Proposition 1 allows us to explicitly express the dependence of the
constant solution $u$ on the parameter $\eta_{0}$ (albeit, in a parametric
form). The same applies to the dependence of $u$ on $\kappa_{\mathrm{s}}$, if
$\hat{W}_{\mathrm{s},0}\neq 0$.

####  IV.1.2 Linear stability

Suppose that $u_{0}\in\mathbb{C}$ is a constant solution of Eq. (9). To
analyze its stability, we insert the ansatz $u(x,t)=u_{0}+v(x,t)$ into Eq. (9)
and linearize the resulting equation with respect to small perturbations
$v(x,t)$. This yields a linear integro-differential equation

$\frac{\displaystyle\partial v}{\displaystyle\partial t}=\mu
v+\frac{\displaystyle\kappa_{\mathrm{v}}}{\displaystyle
2}\mathcal{K}_{\mathrm{v}}(v-\overline{v})+\frac{\displaystyle
i\kappa_{\mathrm{s}}}{\displaystyle
2\pi}\mathcal{K}_{\mathrm{s}}(v+\overline{v}),$ (14)

where

$\mu=-\kappa_{\mathrm{v}}-2iu_{0}.$ (15)

###### Remark 1

If $\mathrm{Re}\>(u_{0})\geq 0$ and hence $u_{0}$ is determined by Proposition
1, then $\mathrm{Re}\>(\mu)<0$. Indeed, this inequality is easy to verify by
inserting the $u_{0}$ defined by formula (10) into (15).

To investigate the growth or decay of different spatial modes, we insert the
ansatz

$v(x,t)=v_{+}(x)e^{\lambda t}+\overline{v}_{-}(x)e^{\overline{\lambda}t}$

into Eq. (14) and equate separately the terms at $e^{\lambda t}$ and
$e^{\overline{\lambda}t}$. Thus, we obtain a spectral problem

$\displaystyle\lambda v_{+}$ $\displaystyle=$ $\displaystyle\mu
v_{+}+\frac{\displaystyle\kappa_{\mathrm{v}}}{\displaystyle
2}\mathcal{K}_{\mathrm{v}}(v_{+}-v_{-})+\frac{\displaystyle
i\kappa_{\mathrm{s}}}{\displaystyle
2\pi}\mathcal{K}_{\mathrm{s}}(v_{+}+v_{-}),$ (16) $\displaystyle\lambda v_{-}$
$\displaystyle=$
$\displaystyle\overline{\mu}v_{-}-\frac{\displaystyle\kappa_{\mathrm{v}}}{\displaystyle
2}\mathcal{K}_{\mathrm{v}}(v_{+}-v_{-})-\frac{\displaystyle
i\kappa_{\mathrm{s}}}{\displaystyle
2\pi}\mathcal{K}_{\mathrm{s}}(v_{+}+v_{-}).$ (17)

Now, we analyze the properties of spatial Fourier modes

$(v_{+}(x),v_{-}(x))^{\mathrm{T}}=(V_{+},V_{-})^{\mathrm{T}}e^{imx}\quad\mbox{with}\quad
m\in\mathbb{Z}\quad\mbox{and}\quad(V_{+},V_{-})^{\mathrm{T}}\in\mathbb{C}^{2}.$

Inserting this ansatz into Eqs. (16), (17) and using the identities

$\mathcal{K}_{\mathrm{v}}e^{imx}=2\pi\hat{W}_{\mathrm{v},m}e^{imx}\quad\mbox{and}\quad\mathcal{K}_{\mathrm{s}}e^{imx}=2\pi\hat{W}_{\mathrm{s},m}e^{imx}$

with the coefficients $\hat{W}_{\mathrm{v},m}$ and $\hat{W}_{\mathrm{s},m}$
defined by (2), we obtain a linear system

$\displaystyle\lambda V_{+}$ $\displaystyle=$ $\displaystyle\mu
V_{+}+\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m}(V_{+}-V_{-})+i\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},m}(V_{+}+V_{-}),$
$\displaystyle\lambda V_{-}$ $\displaystyle=$
$\displaystyle\overline{\mu}V_{-}-\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m}(V_{+}-V_{-})-i\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},m}(V_{+}+V_{-}).$

Its characteristic equation reads

$\lambda^{2}-2\lambda\mathrm{Re}\>(\mu+P_{m})+|\mu+P_{m}|^{2}-|P_{m}|^{2}=0$

where

$P_{m}=\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m}+i\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},m}.$

Solving this equation, we obtain

$\displaystyle\lambda_{\pm}$ $\displaystyle=$
$\displaystyle\mathrm{Re}\>(\mu+P_{m})\pm\sqrt{|P_{m}|^{2}-\left(\mathrm{Im}\>(\mu+P_{m})\right)^{2}}$
$\displaystyle=$
$\displaystyle\mathrm{Re}\>(\mu)+\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m}\pm\sqrt{\pi^{2}\kappa_{\mathrm{v}}^{2}\hat{W}_{\mathrm{v},m}^{2}+\kappa_{\mathrm{s}}^{2}\hat{W}_{\mathrm{s},m}^{2}-\left(\mathrm{Im}\>(\mu)+\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},m}\right)^{2}}.$

Since the Fourier coefficients of $W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$
always satisfy the decay property
$|\hat{W}_{\mathrm{v},m}|,|\hat{W}_{\mathrm{s},m}|\to 0$ for $|m|\to\infty$,
this means $\lambda_{\pm}\to\mu$ or $\lambda_{\pm}\to\overline{\mu}$ for
$|m|\to\infty$. In other words, any potential instability of $u_{0}$ can be
associated with relatively small indices $m$ in the above formula for
$\lambda_{\pm}$.

Instabilities of the uniform state for coupling functions (7), (8) and
$m=0,1,2,3$ are shown in Fig. 7. Their position in the
$(\kappa_{\mathrm{v}},\kappa_{\mathrm{s}})$-plane indicates that only the
modes with $m=0$ and $m=2$ determine stability boundaries. Note that along the
curve $m=0$, the critical eigenvalues everywhere have non-zero imaginary parts
(at least for the range $\kappa_{\mathrm{s}}\in[0,30]$). In contrast, on the
curve $m=2$ we find this property only for $\kappa_{\mathrm{s}}<13.0$. Above
this value, the pair of complex conjugate eigenvalues turns into a pair of
real eigenvalues, so that the corresponding stability boundary has a corner
point. As a result, two qualitatively different bifurcation scenarios of the
uniform state are expected for $\kappa_{\mathrm{s}}<13.0$ and
$\kappa_{\mathrm{s}}>13.0$. Indeed, for $\kappa_{\mathrm{s}}=10$ the uniform
state with $\mathrm{Re}\>(u)\geq 0$ is stable to the left of the $m=0$ curve
and as $\kappa_{\mathrm{v}}$ increases, a Hopf bifurcation occurs, which leads
to the creation of a time-dependent (but still spatially uniform) solution.
However, for $\kappa_{\mathrm{s}}=20$ the uniform state loses stability to the
$m=2$ mode as $\kappa_{\mathrm{v}}$ is increased. The eigenvalues
corresponding to this instability are real, so this is a Turing bifurcation
creating a stationary spatial pattern.

Figure 7:  (a) Stability boundaries of the uniform state of Eq. (3) determined
by the condition $\max(\mathrm{Re}\>\lambda_{+},\mathrm{Re}\>\lambda_{-})=0$
for $m=0,1,2$ and $3$. The uniform state is stable in the shaded region. (b)
Enlargement of the range around $\kappa_{\mathrm{v}}=1$. Note that as
$\kappa_{s}$ increases above $13.0$, the pair of complex conjugate eigenvalues
for $m=2$ turns into a pair of real eigenvalues, so that the corresponding
stability boundary has a corner point. Other parameters: $\eta_{0}=1$ and
$\gamma=0.5$.

###  IV.2 Stationary states

Stationary spatial patterns are non-constant time-independent solutions of Eq.
(3), see Fig. 5(a). We can find them by solving a discretized version of Eq.
(3). Alternatively, we can find them using a self-consistency equation we are
about to derive. For this, we insert the ansatz $u=a(x)$ into Eq. (3) and
obtain

$\gamma+iF(x)-\kappa_{\mathrm{v}}a(x)-ia^{2}(x)=0,$ (18)

where

$F(x)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(a)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(a).$
(19)

Using Proposition 50, we can solve Eq. (18) for each $x$ separately

$a(x)=\mathcal{S}(F(x),\kappa_{\mathrm{v}},\gamma).$ (20)

Then, inserting the result into Eq. (19), we obtain a self-consistency
equation

$F(x)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>\left(\mathcal{S}\left(F(x),\kappa_{\mathrm{v}},\gamma\vphantom{\sum}\right)\right)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>\left(\mathcal{S}\left(F(x),\kappa_{\mathrm{v}},\gamma\vphantom{\sum}\right)\right).$
(21)

Note that Eq. (21) should be solved with respect to the unknown function
$F(x)$, which, in turn, determines the corresponding stationary state $a(x)$
by formula (20). By construction this procedure allows us to find all
stationary states satisfying $|a(x)|\leq 1$, while every unphysical state that
violates this inequality for at least some $x$ is automatically discarded.

In general, Eq. (21) cannot be solved exactly but only approximately. Given
that $F(x)$ is real and $2\pi$-periodic, it is convenient to approximate it
using a truncated Fourier series

$F(x)=\sum\limits_{m=0}^{2M}\hat{f}_{m}\phi_{m}(x)\quad\mbox{with}\quad\hat{f}_{m}\in\mathbb{R},$
(22)

where

$\displaystyle\phi_{0}(x)$ $\displaystyle=1,$ $\displaystyle\phi_{m}(x)$
$\displaystyle=\sqrt{2}\cos((m+1)x/2)\>\>\mbox{for odd}\>\>m,$ (23)
$\displaystyle\phi_{m}(x)$ $\displaystyle=\sqrt{2}\sin(mx/2)\>\>\mbox{for
even}\>\>m,$

are elements of the $L^{2}$-orthogonal trigonometric basis. Inserting this
ansatz into Eq. (21) and projecting the resulting relation onto the basis
function $\phi_{m}(x)$, we obtain a system of equations

$\hat{f}_{m}=\eta_{0}\delta_{m0}+\int_{0}^{2\pi}\left[\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m^{\prime}}\mathrm{Im}\>\left(\mathcal{S}\left(F(x),\kappa_{\mathrm{v}},\gamma\vphantom{\sum}\right)\right)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\hat{W}_{\mathrm{s},m^{\prime}}\mathrm{Re}\>\left(\mathcal{S}\left(F(x),\kappa_{\mathrm{v}},\gamma\vphantom{\sum}\right)\right)\right]\phi_{m}(x)dx$
(24)

for $m=0,1,2,\dots 2M$, where $\delta_{mn}$ is the Kronecker delta,
$\hat{W}_{\mathrm{v},m}$ and $\hat{W}_{\mathrm{s},m}$ are the Fourier
coefficients of $W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$ given by (2), and

$m^{\prime}=0\>\>\mbox{if}\>\>m=0,\quad
m^{\prime}=(m+1)/2\>\>\mbox{if}\>\>m\>\>\mbox{is odd},\quad
m^{\prime}=m/2\>\>\mbox{if}\>\>m\>\>\mbox{is even}.$

For completeness, we list here four identities

$\displaystyle\int_{0}^{2\pi}(\mathcal{K}_{\mathrm{v}}u)(x)\cos(mx)dx=2\pi\hat{W}_{\mathrm{v},m}\int_{0}^{2\pi}u(x)\cos(mx)dx,$
(25)
$\displaystyle\int_{0}^{2\pi}(\mathcal{K}_{\mathrm{v}}u)(x)\sin(mx)dx=2\pi\hat{W}_{\mathrm{v},m}\int_{0}^{2\pi}u(x)\sin(mx)dx,$
(26)
$\displaystyle\int_{0}^{2\pi}(\mathcal{K}_{\mathrm{s}}u)(x)\cos(mx)dx=2\pi\hat{W}_{\mathrm{s},m}\int_{0}^{2\pi}u(x)\cos(mx)dx,$
(27)
$\displaystyle\int_{0}^{2\pi}(\mathcal{K}_{\mathrm{s}}u)(x)\sin(mx)dx=2\pi\hat{W}_{\mathrm{s},m}\int_{0}^{2\pi}u(x)\sin(mx)dx,$
(28)

which are valid for arbitrary real function $u(x)$, any integer $m$, and any
integral operators $\mathcal{K}_{\mathrm{v}}$ and $\mathcal{K}_{\mathrm{s}}$
with symmetric (even) kernels $W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$,
respectively. We used these identities to obtain formula (24).

Note that for symmetric kernels $W_{\mathrm{v}}(x)$ and $W_{\mathrm{s}}(x)$
the dimensionality of system (24) can be halved. Indeed, it is easy to see
that in this case the operators $\mathcal{K}_{\mathrm{v}}$ and
$\mathcal{K}_{\mathrm{s}}$ are invariant on the subspace of even functions.
Therefore, we can assume that $F(x)$ is even and hence $\hat{f}_{m}=0$ for all
even positive $m$. Importantly, after such an assumption we do not need to
care about the translational symmetry of Eq. (21), which would otherwise
require the addition of a pinning condition to (24).

In conclusion, we make a remark that for each stationary pattern $a(x)$ found
using the self-consistency equation (21), its linear stability can be analyzed
according to the scheme of Section IV.1.2. For brevity, we do not perform such
an analysis here, but note that conceptually this will not be very different
from the linear stability analysis of bump states described in [36].

###  IV.3 Travelling wave

In this section, we consider the traveling wave solutions of Eq. (3), which
are uniformly drifting spatial patterns, see Fig. 4(a), with a mathematical
form

$u=a(x+st)\quad\mbox{where}\quad s\neq 0.$

Our primary approach to studying them is to realise that such waves are
stationary in a uniformly travelling coordinate system, travelling at the
speed of the wave. Thus they are steady states of

$\frac{\displaystyle\partial a}{\displaystyle\partial
t}=\gamma-\kappa_{\mathrm{v}}a+i\left[\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(a)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(a)-a^{2}\right]-s\frac{\displaystyle\partial
a}{\displaystyle\partial x}$ (29)

where $s$ is the speed at which they travel [28, 26]. Thus we can find them by
uniformly discretising (29) in space and solving the corresponding large set
of coupled algebraic equations using Newton’s method. The spatial derivative
can be approximated using finite differences. One advantage of this method is
that the stability of the wave can be determined from the eigenvalues of the
linearisation of (29) about a steady state.

An alternative approach to the study of travelling waves in Eq. (3) consists
in deriving a self-consistency equation for them. For this, we note that if
$s>0$ (which will be the default case below) and $a(x)$ is a time-independent
solution of Eq. (29), then $a(x)$ satisfies

$\frac{\displaystyle da}{\displaystyle
dx}=\frac{\displaystyle\gamma}{\displaystyle s}+\frac{\displaystyle
iF(x)}{\displaystyle s}-\frac{\displaystyle\kappa_{\mathrm{v}}}{\displaystyle
s}a-\frac{\displaystyle i}{\displaystyle s}a^{2},$ (30)

where

$F(x)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(a)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(a).$

For Eq. (30) we know (see Proposition 2 in Appendix) that for any $\gamma>0$,
$s>0$, any $\kappa_{\mathrm{v}}\in\mathbb{R}$ and any real $2\pi$-periodic
function $F(x)$, this equation has one and only one $2\pi$-periodic solution,
which lies entirely in the right half-plane
$\mathbb{P}=\\{w\in\mathbb{C}\>:\>\mathrm{Re}\>\>w>0\\}$. Therefore, we can
write

$a(x)=\mathcal{U}(F(x),\kappa_{\mathrm{v}},\gamma,s),$ (31)

where $\mathcal{U}(\cdot)$ denotes the corresponding solution operator. To be
consistent with the above definition of $F(x)$, we must have

$F(x)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>\left(\mathcal{U}\left(F(x),\kappa_{\mathrm{v}},\gamma,s\vphantom{\sum}\right)\right)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>\left(\mathcal{U}\left(F(x),\kappa_{\mathrm{v}},\gamma,s\vphantom{\sum}\right)\right).$
(32)

Thus, in order to find a traveling wave solution of Eq. (3), we can first
solve Eq. (32) for $F(x)$ and then use (31) to calculate the corresponding
wave profile. Note that due to the translational symmetry of Eq. (3) we need
to add an additional pinning condition such as

$\int_{0}^{2\pi}F(x)\sin x\>dx=0$ (33)

to pick up a unique solution of this equation. This condition also allows us
to find the speed of traveling wave $s$.

Let us make a few remarks about the numerical implementation of the method
based on the self-consistency equation (32). Recall that the algorithm for
calculating the solution operator $\mathcal{U}(\cdot)$ is described in Remark
2 in Appendix. Moreover, the oprtimal strategy for solving Eq. (32) is the
application of the Galerkin method. For this, we insert the ansatz (22) into
Eq. (32) and project the resulting relation onto the trigonometric basis
$\phi_{m}(x)$ defined in (23). As a result we obtain a system of $2M+1$
equations

$\hat{f}_{m}=\eta_{0}\delta_{m0}+\int_{0}^{2\pi}\left[\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m^{\prime}}\mathrm{Im}\>\left(\mathcal{U}\left(F(x),\kappa_{\mathrm{v}},\gamma,s\vphantom{\sum}\right)\right)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\hat{W}_{\mathrm{s},m^{\prime}}\mathrm{Re}\>\left(\mathcal{U}\left(F(x),\kappa_{\mathrm{v}},\gamma,s\vphantom{\sum}\right)\right)\right]\phi_{m}(x)dx$
(34)

where

$m^{\prime}=0\>\>\mbox{if}\>\>m=0,\quad
m^{\prime}=(m+1)/2\>\>\mbox{if}\>\>m\>\>\mbox{is odd},\quad
m^{\prime}=m/2\>\>\mbox{if}\>\>m\>\>\mbox{is even}.$

Note that the solutions of Eq. (32) do not necessarily have to be reflection-
symmetric, so all integer indices $m=0,1,\dots,2M$ must be taken into account
in system (34). Moreover, in this case the pinning condition (33) cannot be
omitted, although it becomes particularly simple

$\hat{f}_{2}=0.$ (35)

Combining (34) with (35), we obtain a system that determines all $\hat{f}_{m}$
in (22) and the wave speed $s$.

###  IV.4 Standing wave

We use the term standing wave to refer to a spatiotemporal pattern that is
non-constant in the spatial direction and exhibits periodic oscillation with
the same collective period for each position $x$. Examples of such patterns
are shown in Fig. 1(a), Fig. 3 and Fig. 5(b),(c). Mathematically, each
standing wave corresponds to spatially heterogeneous periodic solution of Eq.
(3), i.e. to a solution of the form $u=a(x,t)$ which is $T$-periodic in $t$
for some $T>0$. To study how these periodic solutions depend on system
parameters, we can again use two approaches. On the one hand, we can look for
periodic solutions of the discretized version of Eq. (3). On the other hand,
we can derive a self-consistency equation similar to Eq. (21) and Eq. (32). To
write such a self-consistency equation, we first perform time-rescaling, by
defining a new unknown function $v(x,t)=u(x,t/\omega)$ where $\omega=2\pi/T$
is the cyclic frequency corresponding to period $T$. The new function $v(x,t)$
satisfies

$\omega\frac{\displaystyle\partial v}{\displaystyle\partial
t}=\gamma-\kappa_{\mathrm{v}}v+i\left[\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(v)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(v)-v^{2}\right],$

or equivalently

$\frac{\displaystyle\partial v}{\displaystyle\partial
t}=\frac{\displaystyle\gamma}{\displaystyle\omega}+\frac{\displaystyle
iF(x,t)}{\displaystyle\omega}-\frac{\displaystyle\kappa_{\mathrm{v}}}{\displaystyle\omega}v-\frac{\displaystyle
i}{\displaystyle\omega}v^{2},$ (36)

where

$F(x,t)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>(v)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>(v).$
(37)

For every fixed $x\in[0,2\pi]$, Eq. (36) has the same form as Eq. (30), but
with other variables $x\mapsto t$ and $s\mapsto\omega$. Therefore, using the
notation of the solution operator $\mathcal{U}(\cdot)$, we can write

$v(x,t)=\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\right).$
(38)

Inserting this into the definition of $F(x,t)$, we obtain a self-consistency
equation

$F(x,t)=\eta_{0}+\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right)+\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right).$
(39)

Recalling that every solution $v_{0}(x,t)$ of the periodic boundary value
problem for Eq. (36) determines a two-parameter familily of solutions
$v_{0}(x+x_{0},t+t_{0})$ with $x_{0},t_{0}\in(0,2\pi)$, we equip Eq. (39) with
two additional pinning conditions

$\int_{0}^{2\pi}\int_{0}^{2\pi}F(x,t)\sin
x\>dx\>dt=\int_{0}^{2\pi}\int_{0}^{2\pi}F(x,t)\sin t\>dx\>dt=0,$ (40)

which allow us to pick up a single function from this family.

In the following we use a Fourier sum approximation for $F(x,t)$:

$F(x,t)=\sum\limits_{m=0}^{2M}\sum\limits_{n=0}^{2N}\hat{f}_{mn}\phi_{m}(x)\phi_{n}(t),$
(41)

where $\hat{f}_{mn}\in\mathbb{R}$, and the $\phi_{m}$ are given in (23).
Applying the standard Galerkin scheme to Eq. (39) we find a system of
algebraic equations that determines the coefficients $\hat{f}_{mn}$:

$\displaystyle\hat{f}_{mn}=\eta_{0}\delta_{m0}\delta_{n0}$ $\displaystyle+$
$\displaystyle\left\langle\kappa_{\mathrm{v}}\mathcal{K}_{\mathrm{v}}\mathrm{Im}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right)\right.$
(42) $\displaystyle+$
$\displaystyle\left.\frac{\displaystyle\kappa_{\mathrm{s}}}{\displaystyle\pi}\mathcal{K}_{\mathrm{s}}\mathrm{Re}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right),\phi_{m}(x)\phi_{n}(t)\right\rangle$

where

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

Moreover, using the identities (25)–(28), we simplify Eq. (42) to obtain

$\displaystyle\hat{f}_{mn}=\eta_{0}\delta_{m0}\delta_{n0}$ $\displaystyle+$
$\displaystyle\left\langle
2\pi\kappa_{\mathrm{v}}\hat{W}_{\mathrm{v},m^{\prime}}\mathrm{Im}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right)\right.$
(43) $\displaystyle+$
$\displaystyle\left.2\kappa_{\mathrm{s}}\hat{W}_{\mathrm{s},m^{\prime}}\mathrm{Re}\>\left(\mathcal{U}\left(F(x,t),\kappa_{\mathrm{v}},\gamma,\omega\vphantom{\sum}\right)\right),\phi_{m}(x)\phi_{n}(t)\right\rangle,$

where

$m^{\prime}=0\>\>\mbox{if}\>\>m=0,\quad
m^{\prime}=(m+1)/2\>\>\mbox{if}\>\>m\>\>\mbox{is odd},\quad
m^{\prime}=m/2\>\>\mbox{if}\>\>m\>\>\mbox{is even}.$

For the patterns shown in Fig. 1(a) and Fig. 3, we can, without loss of
generality, assume $u(-x,t)=u(x,t)$ and hence $F(-x,t)=F(x,t)$. In this case,
$\hat{f}_{mn}=0$ for all even positive $m$. Therefore, the first pinning
condition in (40) is satisfied automatically. For the second pinning
condition, inserting the ansatz (41) into it, we obtain

$\hat{f}_{02}=0.$ (44)

In summary, solving the system (43), (44) we are able to find an approximate
solution $F(x,t)$ of the form (41) as well as the corresponding frequency
$\omega$. Then, by formula (38) we can calculate the solution of Eq. (36),
which after the appropriate time rescaling $u(x,t)=v(x,\omega t)$ yields the
standing wave solution of Eq. (3).

###  IV.5 Lurching solutions

In this section, we consider lurching waves — the most complex type of
solutions of Eq. (3) that we have been able to investigate analytically. A
lurching wave is a travelling wave, but rather than having a constant profile,
its profile periodically oscillates; see Fig. 1(b), Fig. 4(b) and Fig. 5(d).
These are sometimes referred to as modulated travelling waves [39]. Each such
wave can be characterized as a fixed point of a “shift and run” map of the
type described in [41], but with a continuous shift rather than the discrete
one used in [41]. More precisely, we look for lurching waves as fixed points
of the map

$u_{n+1}=P(u_{n};\chi,T)$

where $P(u;\chi,T)$ is defined by taking the function $u(x)$, shifting it in
space by an amount $\chi$ (to the left if the pattern is moving to the right,
recalling that the domain is periodic) and then integrating (3) with this
shifted pattern as the initial condition for an amount of time $T$. A
continuous shift is easily implemented in Fourier space. If

$u(x)=\sum_{n=-M}^{M}a_{n}e^{inx}$

then

$u(x+\chi)=\sum_{n=-M}^{M}\hat{a}_{n}e^{inx}$

where $\hat{a}_{n}=a_{n}e^{in\chi}$. Note that the stability of a lurching
wave can be found from the eigenvalues of the Jacobian of the map $P$
evaluated at the relevant fixed point [41].

It is natural to expect that lurching waves can also be described using a
variant of the self-consistency approach. However, so far we have not been
able to do this. The main problem is that after time rescaling
$v(x,t)=u(x,t/\omega)$ we need to consider Eq. (36) on the functional space

$\tilde{C}_{\chi}=\left\\{v\in
C([0,2\pi]^{2};\mathbb{C})\>:\>v(x+2\pi,t)=v(x,t),\quad
v(x+\chi,t+2\pi)=v(x,t)\right\\}$

and show that for every $F\in\tilde{C}_{\chi}$ equation (36) has a unique
solution $v\in\tilde{C}_{\chi}$ satisfying $|v|\leq 1$. This turns out to be a
nontrivial task due to the nonlocal character of the boundary condition
$v(x+\chi,t+2\pi)=v(x,t)$. So we leave it for future investigation.

##  V Results

Now we show how the analytical methods developed in Section IV can be used to
explain the empirical stability diagrams from Section III. We begin by
reviewing the computational efficiency of the two approaches described in
Section IV. Then, we show the calculated bifurcation diagrams associated with
Figures 2 and 6, which make visible the relationships between various
spatiotemporal patterns which are solutions of Eq. (3), sometimes via unstable
(and therefore invisible) solutions of this equation.

###  V.1 Computational efficiency of the proposed methods

  * •

Stationary states. To find stationary states of (3) one could uniformly
discretise this in space and solve the corresponding large set of coupled
algebraic equations using Newton’s method. Or one could solve the system (24).
Which method is more efficient depends on the number of spatial points used
(call this $\hat{N}$) and $M$, the number of spatial harmonics used to
approximate $F(x)$. If $\hat{N}$ is small and $M$ is large the former method
is more efficient, but if $\hat{N}$ is large and $M$ is small the latter
method is more efficient. The former method has the advantage that the
stability can be readily found from the eigenvalues of the linearisation about
the stationary state.

  * •

Travelling waves. It was found to be much more efficient to use the approach
discussed at the beginning of Sec. IV.3, i.e. equation (29), than to solve the
system of equations (34). The main reasons for this are (i) the large number
of harmonics used to approximate $F$: we used $M=50$, giving 101 unknowns, and
(ii) the fact that to evaluate the operator $\mathcal{U}$ once we need to
numerically integrate (30) four times, as explained in Remark 2. Note that if
the kernels $W_{\mathrm{v}}$ and $W_{\mathrm{s}}$ were described exactly by a
small number of sinusoidal functions, $F$ could be expressed exactly as a
finite Fourier series, with a corresponding small number of unknowns [35, 29].

  * •

Standing waves. A standing wave is a periodic solution of (3) and so could be
studied using conventional methods after spatial discretisation [26]. This
method was found to be much more efficient than the method proposed in Sec.
IV.4, due to the large number of unknowns required for that method,
$(2M+1)(2N+1)$, and the need to integrate (36) four times, as explained above.

###  V.2 Bifurcation diagrams for $\kappa_{\mathrm{s}}=10$

Recall that we found two coexisting stable solutions at
$\kappa_{\mathrm{v}}=1$. We first discuss the standing waves and other related
solutions of Eq. (3). The spatially-uniform state is stable for
$\kappa_{\mathrm{v}}<0.96934$ and it undergoes a supercritical Hopf
bifurcation with $m=0$ at this point. The emerging uniformly oscillating state
is shown in black in Fig. 8 and is stable upon creation but undergoes a
subcritical bifurcation at $\kappa_{\mathrm{v}}=0.96974$, generating an
unstable branch of standing waves, shown in blue in Fig. 8. (This branch of
spatially uniform oscillating states is stable only over a very small range of
$\kappa_{\mathrm{v}}$ values and was not observed in Fig. 2.) At
$\kappa_{\mathrm{v}}=0.94185$ a saddle-node bifurcation occurs creating stable
and unstable branches of standing waves. At $\kappa_{\mathrm{v}}=0.97069$ the
branch of standing waves undergoes a symmetry breaking bifurcation and a
branch of stable standing waves with more complex behavior and less symmetry
(shown in red Fig. 8) appears. The latter loses its stability at
$\kappa_{\mathrm{v}}=1.0091$. This analysis provide an explanation for the
results shown in the left column of Fig. 2: there is a stable standing wave
for smaller values of $\kappa_{\mathrm{v}}$ and one with less symmetry for
larger values.

Figure 8: Bifurcation diagram of standing waves for $\eta_{0}=1$, $\gamma=0.5$
and $\kappa_{\mathrm{s}}=10$. Black: spatially-uniform states. Blue: standing
waves. Red: standing waves with less symmetry. Solid: stable; dashed:
unstable.

We now discuss the travelling waves. The travelling wave is stable at
$\kappa_{\mathrm{v}}=0.96$ and either decreasing or increasing
$\kappa_{\mathrm{v}}$ destabilises the wave in a Hopf bifurcation; see Fig. 9.
The leftmost bifurcation, at $\kappa_{\mathrm{v}}=0.95243$, seems subcritical
and occurs to the right of the saddle-node bifurcation. Following the branch
through the saddle-node bifurcation it terminates in a collision with the
spatially unform state at $\kappa_{\mathrm{v}}=0.9868$. Put another way, the
spatially uniform state undergoes a subcritical Hopf bifurcation with $m=2$ at
this parameter value, leading to the creation of an unstable travelling wave —
see Fig. 7.

Figure 9: Bifurcation diagram of travelling waves for $\eta_{0}=1$,
$\gamma=0.5$ and $\kappa_{\mathrm{s}}=10$. Solid: stable; dashed: unstable.
The two bifurcations at which the wave loses stability are Hopf bifurcations.

The rightmost Hopf bifurcation of the travelling wave, at
$\kappa_{\mathrm{v}}=0.96398$, is supercritical and results in the creation of
a stable lurching solution. We followed this lurching wave and a plot of the
spatial shift $\chi$ versus $\kappa_{\mathrm{v}}$ is shown in Fig. 10. The
magnitude of $\chi$ is small, in keeping with the slow drift mentioned in Sec.
III.1. These lurching waves remain stable until $\kappa_{\mathrm{v}}=1.0301$,
where they apparently undergo a torus bifurcation, creating a solution which
is quasiperiodic in a uniformly travelling coordinate frame. As
$\kappa_{\mathrm{v}}$ grows further, these waves become chaotic. This analysis
provide an explanation for the results shown in the right column of Fig. 2:
the travelling wave in Fig. 9 has $\Delta R=0$, the lurching wave is seen at
$\kappa_{\mathrm{v}}=1$, for example, and the quasiperiodic and chaotic
behaviour is seen in Fig. 2(e).

Figure 10: Behaviour of lurching waves for $\eta_{0}=1$, $\gamma=0.5$ and
$\kappa_{\mathrm{s}}=10$. $\chi$ is the spatial shift. Solid: stable; dashed:
unstable. The left end of the branch of solutions is created in the rightmost
bifurcation shown in Fig. 9.

###  V.3 Bifurcation diagram for $\kappa_{\mathrm{s}}=20$.

The uniform state undergoes a subcritical Turing bifurcation with $m=2$ at
$\kappa_{\mathrm{v}}=-1.53$ (see Fig. 7). An unstable branch of stationary
two-bump solutions emerges at this point, shown dashed in Fig. 11. This
solution becomes stable in a saddle-node bifurcation at
$\kappa_{\mathrm{v}}=-1.6099$. At $\kappa_{\mathrm{v}}=0.88565$ a
supercritical Hopf bifurcation occurs and a branch of standing waves (or
breathing two-bump solutions), shown in blue in Fig. 12 emerges.

Figure 11: Stationary two-bump solution $\eta_{0}=1$, $\gamma=0.5$ and
$\kappa_{\mathrm{s}}=20$. Solid: stable; dashed: unstable. The vertical axis
shows $\max_{x}R(x)-\min_{x}R(x)$. The solution is stable between the saddle-
node bifurcation and a supercritical Hopf bifurcation at
$\kappa_{\mathrm{v}}=0.88565$. (Note that negative values of
$\kappa_{\mathrm{v}}$ are unphysical, but we plot them to show the full branch
of solutions.)

These standing waves are stable until a symmetry breaking period-doubling
bifurcation at $\kappa_{\mathrm{v}}=1.0719$. To the right from this point, a
branch of standing waves in the form of alternating two-bump patterns is
stable, shown in solid red in Fig. 12. At $\kappa_{\mathrm{v}}=1.5563$ this
alternating pattern becomes unstable through a supercritical torus bifurcation
and a branch of lurching states appears, shown in black in Fig. 12. The
lurching state is stable until $\kappa_{\mathrm{v}}=1.6248$. After that it
undergoes a sequence of bifurcations leading to a chaotic regime. Note that
the branch of periodic solutions shown in blue in Fig. 12 terminates in a
collision with an unstable spatially uniform periodic solution (not shown or
analysed). This analysis provides an explanation for the results shown in Fig.
6.

Figure 12: Bifurcation diagrams for $\kappa_{\mathrm{s}}=20$. Blue: periodic
solution. Red: period-doubled solution. Black: lurching solution. Solid:
stable; dashed: unstable. Panel (b) is an enlargment of panel (a). Panel (c)
shows the shift, $\chi$, of the lurching solution.

##  VI Discussion

In this paper, we considered the neural field equation (3), which describes
the long-term dynamics of a ring network of QIF neurons with nonlocal synaptic
and gap junction coupling, in the limit of a large number of neurons. We
showed that apart from spatially uniform states, this model also supports
various spatiotemporal patterns, including stationary bump states, standing
and travelling waves, and lurching waves. We showed how each of these states
can be analyzed semi-analytically using the discretized version of Eq. (3) or
using the corresponding self-consistency equation. We believe this to be the
first study of QIF networks using self-consistency arguments. Finally, we
computed detailed bifurcation diagrams for (3), which linked the previous
observations together.

The linear stability analysis of constant solutions of Eq. (3) performed in
Sec. IV.1 revealed that there are two qualitatively different bifurcation
scenarios for small and large values of synaptic coupling; see Fig. 7. More
specifically, for $\kappa_{\mathrm{s}}<13$ the uniform state loses its
stability due to the Hopf bifurcation with a spatially uniform eigenfunction,
resulting in the creation of a spatially-uniform oscillating state. In
contrast, for $\kappa_{\mathrm{s}}>13$ the destabilization of the uniform
state occurs due to a Turing bifurcation with a spatially modulated
eigenfunction with wave number $m=2$. This results in the creation of a “two-
bump” stationary solution.

Using the algorithms described in Sec. IV, we followed the five types of
solutions of Eq. (3) primarily as the parameter $\kappa_{\mathrm{v}}$ was
varied. This allowed us to better understand the relationship between these
solutions and provide an explanation of the numerical results found in Sec.
III. More specifically, for $\kappa_{\mathrm{s}}=10$ we found coexisting
standing waves and travelling waves, and a type of lurching wave that resulted
from a travelling wave undergoing an oscillatory instability. For
$\kappa_{\mathrm{s}}=20$ we found the scenario shown in Fig. 12, where more
complex solutions are created as $\kappa_{\mathrm{v}}$ is increased.
Interestingly, the lurching solution found for this value of
$\kappa_{\mathrm{s}}$ results from a standing wave solution losing stability
and starting to travel, hence the increase in $\chi$ from zero as
$\kappa_{\mathrm{v}}$ is increased, shown in Fig. 12(c).

The model (3) has many parameters and we have only investigated varying one of
them for several values of another. However, our methods could equally-well be
used to investigate other variations of parameters. We believe that the
numerical and analytical methods of Sec. III and Sec. IV can be extended to
other types of networks consisting of QIF neurons, for which stationary,
oscillating and moving patterns also define their typical dynamical behaviour.
In particular, these methods can be useful for studying two-, three- and
higher-dimensional arrays, networks with population structure and networks
with propagation delays. We leave these questions for future work.

## Appendix

Let us consider a complex Riccati equation

$\frac{\displaystyle dw}{\displaystyle dt}=\gamma+if(t)-g(t)w-iw^{2}$ (45)

with a parameter $\gamma\in\mathbb{R}$ and real-valued functions $f(t)$ and
$g(t)$. In this section, we show that for $\gamma>0$ all solutions of Eq. (45)
have special invariance properties with respect to the complex half-plane

$\mathbb{P}=\\{w\in\mathbb{C}\>:\>\mathrm{Re}\>\>w>0\\}.$

In order to proceed, we recall that the complex plane transformation

$z=\frac{\displaystyle 1-w}{\displaystyle 1+w}$ (46)

maps the half-plane $\mathbb{P}$ onto the unit disc

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

Moreover, using formula (46) as a change of variable, we can rewrite Eq. (45)
in an equivalent form

$\frac{\displaystyle dz}{\displaystyle
dt}=-\frac{\displaystyle\gamma+if(t)}{\displaystyle
2}(1+z)^{2}+\frac{\displaystyle g(t)}{\displaystyle
2}(1-z^{2})+\frac{\displaystyle i}{\displaystyle 2}(1-z)^{2}.$ (47)

Then, it is easy to see that Eq. (47) is a complex Riccati equation

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

with

$c_{0}(t)=\frac{\displaystyle-\gamma-if(t)+g(t)+i}{\displaystyle 2},\quad
c_{1}(t)=-\gamma-if(t)-i,\quad c_{2}(t)=\frac{\displaystyle-\gamma-
if(t)-g(t)+i}{\displaystyle 2}.$

Therefore, Proposition 2 from [29] implies that for every $\gamma>0$ and for
all real $2\pi$-periodic continuous functions $f(t)$ and $g(t)$, Eq. (47) has
exactly one stable $2\pi$-periodic solution and this solution satisfies
$|z(t)|<1$ for all $0\leq t\leq 2\pi$. Moreover, using the main argument of
the proof of Proposition 2 in [29], we can show that every solution of Eq.
(47) with an initial condition in $\mathbb{D}$ remains in this unit disc for
any finite time interval.

Due to the equivalence between the solutions of Eq. (45) and the solutions of
Eq. (47), we can reformulate the above results as follows.

###### Proposition 2

Suppose that $\gamma>0$ and $f(t)$, $g(t)$ are real $2\pi$-periodic continuous
functions.

(i) Then for every $w_{0}\in\mathbb{P}$ the complex Riccati equation (45) with
initial condition $w(0)=w_{0}$ has a solution $w(t)$, which is bounded on the
interval $[0,2\pi]$ and satisfies $\mathrm{Re}\>\>w(t)>0$.

(ii) Moreover, Eq. (45) has exactly one stable $2\pi$-periodic solution and
this solution satisfies $\mathrm{Re}\>\>w(t)>0$ for all $0\leq t\leq 2\pi$.

Since Eq. (45) belongs to the class of complex Riccati equations its Poincaré
map coincides with a Möbius transformation, see [8, 42]. This fact can be used
to speed up the computation of periodic solutions of Eq. (45).

###### Remark 2

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

(i) One solves Eq. (45) on the interval $t\in(0,2\pi]$ with three different
initial conditions $w(0)=w_{k}\in\mathbb{P}$, $k=1,2,3$, and obtains three
solutions $W_{k}(t)$. Since each $w_{k}$ lies in the half-plane $\mathbb{P}$
this automatically implies that $W_{k}(t)$ is bounded on the interval
$t\in[0,2\pi]$. Note that this fact is not obvious a priori, since some of the
solutions of Eq. (45) blow up in finite time.

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

$\mathcal{M}(w)=\frac{aw+b}{cw+d}$

where

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

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

$cw^{2}+dw-aw-b=0.$

This yields two roots

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

(iv) Choosing from the roots $w_{+}$ and $w_{-}$ the one that lies in the
half-plane $\mathbb{P}$, one obtains the initial condition that determines the
periodic solution of interest. The latter can be computed by solving Eq. (45)
with this initial condition.

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

$|\zeta_{1}-\zeta_{2}|+|\zeta_{3}-\zeta_{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 $(\zeta_{1}+\zeta_{2}+\zeta_{3})/3$.

The above steps (i)–(v) can be understood as a constructive definition of the
solution operator of Eq. (45), which for every $\gamma>0$ and all real
$2\pi$-periodic coefficients $f(t)$ and $g(t)$ yields the corresponding stable
$2\pi$-periodic solution of Eq. (45).

###### Proposition 3

For constant coefficients $f$ and $g$, Eq. (45) has a single fixed point
satisfying $\mathrm{Re}\>(w)\geq 0$. This fixed point is determined by the
formula

$w=\mathcal{S}(f,g,\gamma):=\frac{\displaystyle 1}{\displaystyle
2}\sqrt{\frac{\displaystyle\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}-g^{2}+4f}{\displaystyle
2}}+\frac{\displaystyle i}{\displaystyle
2}\left(g-\sqrt{\frac{\displaystyle\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}+g^{2}-4f}{\displaystyle
2}}\right).$ (50)

Proof: Any fixed point of Eq. (45) satisfies

$\gamma+if-gw-iw^{2}=0.$

This quadratic equation has two solutions

$w=-\frac{\displaystyle 1}{\displaystyle
2i}\left(g\pm\sqrt{g^{2}+4i(\gamma+if)}\right)=\frac{\displaystyle
i}{\displaystyle 2}\left(g\pm\sqrt{g^{2}-4f+4i\gamma}\right).$ (51)

Since $\gamma>0$, de Moivre’s formula yields

$g^{2}-4f+4i\gamma=\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}\left(\cos\left(\mathrm{arccot}\frac{\displaystyle
g^{2}-4f}{\displaystyle
4\gamma}\right)+i\sin\left(\mathrm{arccot}\frac{\displaystyle
g^{2}-4f}{\displaystyle 4\gamma}\right)\right).$

Hence

$\sqrt{g^{2}-4f+4i\gamma}=\sqrt[4]{(g^{2}-4f)^{2}+16\gamma^{2}}\left(\cos\left(\frac{\displaystyle
1}{\displaystyle 2}\mathrm{arccot}\frac{\displaystyle g^{2}-4f}{\displaystyle
4\gamma}\right)+i\sin\left(\frac{\displaystyle 1}{\displaystyle
2}\mathrm{arccot}\frac{\displaystyle g^{2}-4f}{\displaystyle
4\gamma}\right)\right).$

Using the half-angle formulas and the identity

$\cos\left(\mathrm{arccot}\frac{\displaystyle g^{2}-4f}{\displaystyle
4\gamma}\right)=\frac{\displaystyle
g^{2}-4f}{\displaystyle\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}},$

we rewrite this as follows

$\sqrt{g^{2}-4f+4i\gamma}=\sqrt{\frac{\displaystyle\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}+g^{2}-4f}{\displaystyle
2}}+i\sqrt{\frac{\displaystyle\sqrt{(g^{2}-4f)^{2}+16\gamma^{2}}-g^{2}+4f}{\displaystyle
2}}.$

Inserting this expression into (51) and choosing the sign of the square root
that ensures the inequality $\mathrm{Re}\>(w)\geq 0$, we obtain formula (50).

## References

  * [1] Shun-ichi Amari.  Dynamics of pattern formation in lateral-inhibition type neural fields.  Biological cybernetics, 27(2):77–87, 1977.
  * [2] M. Bataille-Gonzalez, M. G. Clerc, E. Knobloch, and O. E. Omel’chenko.  Traveling spiral wave chimeras in coupled oscillator systems: emergence, dynamics, and transitions.  New Journal of Physics, 25(10):103023, 2023.
  * [3] Paul C Bressloff.  Spatiotemporal dynamics of continuum neural fields.  Journal of Physics A: Mathematical and Theoretical, 45(3):033001, 2012.
  * [4] Paul C Bressloff and Samuel R Carroll.  Laminar neural field model of laterally propagating waves of orientation selectivity.  PLoS computational biology, 11(10):e1004545, 2015.
  * [5] Paul C Bressloff and Matthew A Webber.  Neural field model of binocular rivalry waves.  Journal of computational neuroscience, 32(2):233–252, 2012.
  * [6] Áine Byrne, Daniele Avitabile, and Stephen Coombes.  Next-generation neural field model: The evolution of synchrony within patterns and waves.  Phys. Rev. E, 99:012313, Jan 2019.
  * [7] Áine Byrne, James Ross, Rachel Nicks, and Stephen Coombes.  Mean-field models for EEG/MEG: from oscillations to waves.  Brain topography, 35(1):36–53, 2022.
  * [8] J. Campos.  Möbius transformations and periodic solutions of complex Riccati equations.  Bull. London Math. Soc., 29:205–215, 1997.
  * [9] Albert Compte, Nicolas Brunel, Patricia S. Goldman-Rakic, and Xiao-Jing Wang.  Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model.  Cerebral Cortex, 10(9):910–923, 2000.
  * [10] S Coombes.  Dynamics of synaptically coupled integrate-and-fire-or-burst neurons.  Physical Review E, 67(4):041910, 2003.
  * [11] S. Coombes.  Waves, bumps, and patterns in neural field theories.  Biol. Cybern., 93(2):91–108, 2005.
  * [12] Stephen Coombes, Peter beim Graben, Roland Potthast, and James Wright, editors.  Neural Fields: Theory and Applications.  Springer, 2014.
  * [13] Stephen Coombes and Kyle CA Wedgwood.  Neurodynamics: An Applied Mathematics Perspective, volume 75.  Springer Nature, 2023.
  * [14] B. Ermentrout.  Neural networks as spatio-temporal pattern-forming systems.  Rep. Prog. Phys., 61:353–430, 1998.
  * [15] Bard Ermentrout, Jonathan Rubin, and Remus Osan.  Regular traveling waves in a one-dimensional network of theta neurons.  SIAM Journal on Applied Mathematics, 62(4):1197–1221, 2002.
  * [16] G.B. Ermentrout and D. Terman.  Mathematical Foundations of Neuroscience, volume 35.  Springer Verlag, 2010.
  * [17] Jose M. Esnaola-Acebes, Alex Roxin, Daniele Avitabile, and Ernest Montbrió.  Synchrony-induced modes of oscillation of a neural field model.  Phys. Rev. E, 96:052407, Nov 2017.
  * [18] Stefanos E Folias and Paul C Bressloff.  Breathing pulses in an excitatory neural network.  SIAM Journal on Applied Dynamical Systems, 3(3):378–407, 2004.
  * [19] David Golomb and G Bard Ermentrout.  Continuous and lurching traveling pulses in neuronal networks with delay and spatially decaying connectivity.  Proceedings of the National Academy of Sciences, 96(23):13480–13485, 1999.
  * [20] David Golomb and G Bard Ermentrout.  Effects of delay on the type and velocity of travelling pulses in neuronal networks with spatially decaying connectivity.  Network: Computation in Neural Systems, 11(3):221, 2000.
  * [21] David Golomb and G Bard Ermentrout.  Bistability in pulse propagation in networks of excitatory and inhibitory populations.  Physical Review Letters, 86(18):4179, 2001.
  * [22] David Horn and Irit Opher.  Solitary waves of integrate-and-fire neural fields.  Neural Comput., 9(8):1677–1690, 1997.
  * [23] Xiaoying Huang, William C Troy, Qian Yang, Hongtao Ma, Carlo R Laing, Steven J Schiff, and Jian-Young Wu.  Spiral waves in disinhibited mammalian neocortex.  Journal of Neuroscience, 24(44):9897–9902, 2004.
  * [24] Mikail Khona and Ila R Fiete.  Attractor and integrator networks in the brain.  Nature Reviews Neuroscience, 23(12):744–766, 2022.
  * [25] Carlo R Laing.  Derivation of a neural field model from a network of theta neurons.  Physical Review E, 90(1):010901, 2014.
  * [26] Carlo R Laing.  Numerical bifurcation theory for high-dimensional neural models.  The Journal of Mathematical Neuroscience, 4(1):13, 2014.
  * [27] Carlo R Laing.  Exact neural fields incorporating gap junctions.  SIAM Journal on Applied Dynamical Systems, 14(4):1899–1929, 2015\.
  * [28] Carlo R Laing and Oleh Omel’chenko.  Moving bumps in theta neuron networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(4):043117, 2020.
  * [29] Carlo R Laing and Oleh E Omel’chenko.  Periodic solutions in next generation neural field models.  Biological Cybernetics, 117(4):259–274, 2023.
  * [30] Carlo R Laing and William C Troy.  Two-bump solutions of amari-type models of neuronal pattern formation.  Physica D: Nonlinear Phenomena, 178(3-4):190–218, 2003.
  * [31] Carlo R Laing, William C Troy, Boris Gutkin, and G Bard Ermentrout.  Multiple bumps in a neuronal model of working memory.  SIAM Journal on Applied Mathematics, 63(1):62–97, 2002.
  * [32] Ernest Montbrió and Diego Pazó.  Exact mean-field theory explains the dual role of electrical synapses in collective synchronization.  Physical review letters, 125(24):248101, 2020.
  * [33] Ernest Montbrió, Diego Pazó, and Alex Roxin.  Macroscopic description for networks of spiking neurons.  Phys. Rev. X, 5:021028, Jun 2015.
  * [34] O. E. Omel’chenko.  Mathematical framework for breathing chimera states.  Journal of Nonlinear Science, 32(2):1–34, 2022.
  * [35] O. E. Omel’chenko.  Periodic orbits in the ott-antonsen manifold.  Nonlinearity, 36:845–861, 2023.
  * [36] Oleh Omel’chenko and Carlo R Laing.  Collective states in a ring network of theta neurons.  Proceedings of the Royal Society A, 478(2259):20210817, 2022.
  * [37] Bastian Pietras, Federico Devalle, Alex Roxin, Andreas Daffertshofer, and Ernest Montbrió.  Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks.  Physical Review E, 100(4):042412, 2019.
  * [38] David J Pinto and G Bard Ermentrout.  Spatially structured activity in synaptically coupled neuronal networks: I. Traveling fronts and pulses.  SIAM journal on Applied Mathematics, 62(1):206–225, 2001.
  * [39] Björn Sandstede and Arnd Scheel.  On the structure of spectra of modulated travelling waves.  Mathematische Nachrichten, 232(1):39–93, 2001.
  * [40] Helmut Schmidt and Daniele Avitabile.  Bumps and oscillons in networks of spiking neurons.  Chaos, 30:033133, 2020.
  * [41] Thomas M Wasylenko, Jaime E Cisternas, Carlo R Laing, and Ioannis G Kevrekidis.  Bifurcations of lurching waves in a thalamic neuronal network.  Biological cybernetics, 103(6):447–462, 2010.
  * [42] P. Wilczyński.  Planar nonautonomous polynomial equations: the Riccati equation.  J. Differ. Equ., 244:1304–1328, 2008.
  * [43] Hugh R Wilson and Jack D Cowan.  A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue.  Kybernetik, 13(2):55–80, 1973.
  * [44] Klaus Wimmer, Duane Q Nykamp, Christos Constantinidis, and Albert Compte.  Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory.  Nature neuroscience, 17(3):431–439, 2014.
  * [45] Xiaohui Xie, Richard HR Hahnloser, and H Sebastian Seung.  Double-ring network model of the head-direction system.  Physical Review E, 66(4):041902, 2002.
  * [46] Alice C Yew, DH Terman, and G Bard Ermentrout.  Propagating activity patterns in thalamic neuronal networks.  SIAM Journal on Applied Mathematics, 61(5):1578–1604, 2001.