Bumps and waves in spiking networksD. Avitabile, J. L. Davis, K. C. A.
Wedgwood hypothesisHypothesis problemProblem remarkRemark

# Bump attractors and waves in networks of leaky integrate-and-fire neurons

Daniele Avitabile Department of Mathematics, Vrije Universiteit Amsterdam and
MathNeuro Team, Inria Sophia Antipolis. , https://www.danieleavitabile.com.
<EMAIL_ADDRESS>Joshua L. Davis  Defence Science and Technology Laboratory,
Cyber and Information Systems Division. Kyle C. A. Wedgwood  Living Systems
Institute, College of Engineering, Mathematics and Physical Sciences.

###### Abstract

Bump attractors are wandering localised patterns observed in in vivo
experiments of spatially-extended neurobiological networks. They are important
for the brain’s navigational system and specific memory tasks. A bump
attractor is characterised by a core in which neurons fire frequently, while
those away from the core do not fire. These structures have been found in
simulations of spiking neural networks, but we do not yet have a mathematical
understanding of their existence, because a rigorous analysis of the nonsmooth
networks that support them is challenging. We uncover a relationship between
bump attractors and travelling waves in a classical network of excitable,
leaky integrate-and-fire neurons. This relationship bears strong similarities
to the one between complex spatiotemporal patterns and waves at the onset of
pipe turbulence. Waves in the spiking network are determined by a firing set,
that is, the collection of times at which neurons reach a threshold and fire
as the wave propagates. We define and study analytical properties of the
voltage mapping, an operator transforming a solution’s firing set into its
spatiotemporal profile. This operator allows us to construct localised
travelling waves with an arbitrary number of spikes at the core, and to study
their linear stability. A homogeneous “laminar” state exists in the network,
and it is linearly stable for all values of the principal control parameter.
Sufficiently wide disturbances to the homogeneous state elicit the bump
attractor. We show that one can construct waves with a seemingly arbitrary
number of spikes at the core; the higher the number of spikes, the slower the
wave, and the more its profile resembles a stationary bump. As in the fluid-
dynamical analogy, such waves coexist with the homogeneous state, are
unstable, and the solution branches to which they belong are disconnected from
the laminar state; we provide evidence that the dynamics of the bump attractor
displays echoes of the unstable waves, which form its building blocks.

##  1 Introduction

Figure 1: (a) Bifurcation diagram of travelling waves in a continuous
integrate-and-fire model. (b) Bump attractor in a discrete integrate-and-fire
model with 5000 neurons (dots represent neuronal firing events). Model
descriptions and parameters will be given later in Sections 2, 1, 2 and 3).
The bifurcation diagram in (a) shows selected branches of stable (blue) and
unstable (grey) travelling waves, in the continuation parameter $\beta$, that
is, the timescale at which neuron process incoming currents. Waves are
measured using their width $\Delta$ and are indexed by the number of advected
spikes. The profile of $\textrm{TW}_{11}$, a representative wave with $11$
spikes, is shown. A large number of waves
($\textrm{TW}_{2}$–$\textrm{TW}_{160}$ in the picture, but many more unstable
branches are omitted) coexist with the trivial homogeneous state, which is the
only steady state in the model, and which is stable for all values of $\beta$.
Thin waves are stable, with a small basin of attraction. Sufficiently large,
localised disturbances of the homogenous state lead to the formation of a bump
with a characteristic width: the bump in (b) is marked as (2) in (a). The
region in parameter space where bumps are observed is crowded with unstable
travelling waves, with a large number of spikes, and a width comparable to the
one of the bump. Branches of waves are detached from the homogeneous state;
they originate at critical points called grazing points (blue dots in (a));
waves that are born stable become unstable at oscillatory bifurcations (grey
dots in (a)).

Understanding how networks of coupled, excitable units generate collective
patterns is a central question in the life sciences and, more generally, in
applied mathematics. In particular, the study of network models is ingrained
in neuroscience applications, as they provide a natural way to describe the
interaction of neurons within a population, or of neural populations within
the cortex. In the past decades, a large body of work in mathematical
neuroscience has addressed the development and analysis of neurobiological
networks, with the view of studying the origin of large-scale brain activity
[35, 13, 21], and mapping single-cell and population parameters to
experimental observations, including in vivo and in vitro cortical waves [70,
47, 43], electroencephalogram recordings [77], and patterns in the visual
cortex [17].

This paper presents a novel mathematical characterisation of a prominent
example of spatiotemporal pattern in neuroscience applications, and draws an
analogy inspired by recent progress in the fluid-dynamics literature on
transition to turbulence in a pipe [7]. We focus on the so-called bump
attractor111In the neuroscience literature the term bump attractor refers
sometimes to a network producing a localised pattern, as opposed to the
pattern itself. Similarly, some authors use ring attractor for a network with
ring topology, generating a localised activity bump. Here, we use these terms
to refer to patterns, following the standard convention in the dynamical
system literature., a localised pattern of neural activity observable in
experiments and numerical simulations of spatially-extended, neurobiological
networks [68, 87]. Bump attractors have been associated to working memory, the
temporary storage of information in the brain, and experimental evidence
supporting their existence has been found in the navigational systems of rats
[52] and flies [51, 79], and in oculomotor responses of monkeys [85].

In a bump attractor, the neural activity is localised around a particular
position in the network (see Fig. 1(b)) which may encode, for instance, the
animal’s head position. Bumps are elicited by transient localised stimuli,
such as visual cues at a specific locations, but are sustained autonomously by
the network once the stimulus is removed (the network dynamics is _attracted_
to the bump). These coherent structures display a characteristic wandering
motion, and may exhibit discontinuous jumps if the impinging stimulus
undergoes sudden spatial shifts [51].

###  1.1 Model descriptions

Mathematical neuroscience has a long-standing fascination with localised bumps
of activity. Neural field models, which represent the cortex as a continuum,
were introduced in the 1970s, and spatially-localised solutions to these
models appeared already in seminal papers on the subject, by Wilson and Cowan
[84], and by Amari [1]. Since then, many authors have studied localised
solutions in neural fields, addressed the derivation of neural field equations
from first principles, their relevance to a wide variety of neural phenomena,
and their rigorous mathematical treatment. We refer the reader to [35, 13, 21]
for exhaustive introductions on this topic.

Neural fields are integro-differntial equations which model the cortex as an
excitable, spatially-extended medium. Mathematical mechanisms for pattern
formation in neural fields are similar to the ones found in other nonlinear
media, such as reaction-diffusion systems, albeit their analysis requires some
modifications because these models contain nonlocal operators. Stationary
bumps form via instabilities of the homogeneous steady state, and their
profile depends strongly on the coupling, which typically involves excitation
on short spatial scales, and inhibition on longer scales [35, 13, 21]. Neural
fields support travelling bump solutions, as well as wandering bumps. The
latter are obtained in neural fields that incorporate stochastic terms
deriving, for instance, from noisy currents [50, 59].

Neural fields are heuristic, coarse grained models, hence they bypass
microscopic details that are important in bump attractors. For instance, the
neural firing rate, which is an emergent neural property and an observable in
the bump attractor experiments, is a prescribed feature in neural fields,
hardwired in the model through an ad-hoc firing-rate function. On the other
hand, numerical simulations of large networks of Hodgkin–Huxley-type neurons
with realistic biological details can display emergent neural firing, but
their mathematical treatment is challenging, and still under development [4].

_Spiking neural networks_ are intermediate, bottom-up models which couple
neurons with idealised dynamics. The salient feature of spiking models is that
the firing of a neuron is described as an event, and no attempt is made to
model the temporal evolution of the membrane potential during and after the
spike [49, 39, 13]. Spiking neural networks are specified by 3 main
ingredients: (i) an ordinary differential equation (ODE) for the membrane
potential of each neuron; (ii) rules to define the occurrence and effects of a
spike; (iii) the network coupling.

Since the introduction of the first single-cell spiking model by Lapicque
[57], the so-called _leaky integrate-and-fire model_ , more realistic variants
have been proposed, and spiking neural networks have become a widely adopted
tool in theoretical neuroscience [78, 15, 39]. In specific spiking models,
analytical progress has been made for single neurons and spatially-independent
networks using coordinate transformations [34, 62], dimension reduction [58,
63], and probabilistic methods [26] (see also the reviews [72, 9]). Exact
mean-field reductions, amenable to standard pattern-formation analysis, have
been derived in selected spatially-extended networks [53, 36, 16, 74], but
generally the study of bumps in spiking models has been possible only with
numerical simulations [54, 18].

The present paper investigates localised patterns supported in discrete and
continuous networks of globally coupled leaky integrate-and-fire neurons. In
direct numerical simulations, we use a well-known discrete model, proposed by
Laing and Chow [54], whose details will be given later. For now it will
suffice to consider a cursory formulation of the model, simulated in Fig.
1(b). The network describes the idealised, dimensionless voltage dynamics of
$n$ all-to-all coupled neurons, evenly-spaced in a cortex with ring geometry,

(1) $\dot{v}_{i}=-v_{i}+I_{i}(t)+\sum_{j=1}^{n}S_{ij}(v_{j},\beta),\qquad
i=1,\ldots,n.$

The dynamics of the $i$th neuron’s membrane voltage is specified in terms of
an Ohmic leakage current $-v_{i}$, an external current $I_{i}(t)$, and
voltage-dependent currents, received from other neurons via synaptic
connections; the latter currents, indicated by $S_{ij}$, have a characteristic
time scale $\beta$, and are caused by $v_{j}$ crossing a fixed threshold (when
the $j$th neuron fires). After a firing event, marked with a dot in Fig. 1(b),
the neuron’s voltage is instantly reset to a default value, from which it can
evolve again, following an ODE of type Eq. 1. Discrete and continuous networks
of this type are canonical models of neural activity, widely adopted in the
mathematical neuroscience literature [82, 81, 30, 31, 14, 54, 32, 65, 19, 64,
62, 40]. It is now established that such networks supports bump attractors and
localised waves, but an explanation of the mathematical origins of the former
is still lacking.

This paper presents a new approach to the problem, and uncovers a novel
bifurcation structure for localised travelling waves of the network, shedding
light onto the nature of the bump attractor. Our findings suggest an
intriguing analogy between the bump attractor in the integrate-and-fire
network and the phenomenon of transition to turbulence in a pipe. The analogy
between the bifurcation scenarios of these two problems is notable, and we use
it here to summarise our results, highlighting similarities between the
respective bifurcation structures and dynamical regimes.

###  1.2 Transition to turbulence in a pipe

Stemming from the pioneering experiments of Reynolds [69], a large body of
work in fluid dynamics has addressed how high-speed pipe flows transition from
a laminar state, whose analytical expression is known in closed form, to
complex spatio-temporal patterns, characteristic of the turbulent regime (see
[7] for a recent review). In this context, the Navier-Stokes equations are
studied as a deterministic dynamical system, subject to changes in Reynolds
number, the principal control parameter. Experiments and computer simulations
indicate that the laminar state is stable to infinitesimal perturbations
(linearly stable) up to large values of the control parameter (up to at least
Reynolds number $10^{7}$ in numerical computations) [73, 25, 61, 80, 60].
However, when a disturbance is applied at sufficiently large Reynolds numbers,
a transition to turbulence is observed, depending sensitively on the applied
stimulus [25, 45]. Current opinions view the transition as being determined by
_travelling wave solutions_ to the Navier-Stokes equations [75, 29, 37, 83,
66, 41]. These invariant states, whose spatial profiles display hallmarks of
the turbulent transition, coexist with the laminar state at intermediate
Reynolds numbers, are linearly unstable, and provide an intricate blueprint
for the dynamics, in that orbits may visit transiently these repelling
solutions in phase space. Importantly, the waves lie on branches that are
disconnected from the stable laminar state, and emerge at saddle-node
bifurcations [37, 83]: this turbulence mechanism is therefore different from
other paradigmatic routes to chaos, involving the destabilisation of the
laminar state, and the progressive appearance of more complicated structures
via a cascade of instabilities [56, 46, 71].

###  1.3 Summary of results

In a series of recent papers addressing turbulence from a dynamical-system
viewpoint, Barkley proposed an analogy between pipe flows and excitable media,
using the propagation of an electrical pulse along the axon of a neuron as a
metaphor for localised turbulence puffs [5, 6, 8, 7]. The present paper offers
a specular view, at a different scale: we are motivated by studying a
canonical, complex neurobiological network of coupled excitable neurons,
supporting localised spatio-temporal chaos, and we find a compelling
similarity between the bifurcation structure of waves in this system, and the
one of waves in the pipe turbulence.

With reference to Fig. 1, the principal control parameter of the problem is
$\beta$, the timescale of synaptic currents: a low $\beta$ gives small,
persisting currents, while $\beta\to\infty$ gives large instantaneous
currents. A homogeneous steady state exists and is linearly stable for all
values of $\beta$ ($\Delta=0$ line in Fig. 1(b)), but transient localised
stimuli trigger the bump attractor [54]. In the analogy, the homogeneous
equilibrium plays the role of a “laminar state”. We stress that the
homogeneous steady state is the only equilibrium of the model. Thus, the model
can not support branches of stationary bump solutions. Instead, we demonstrate
that travelling waves are key to understand the bump attractor.

We consider a spatially-continuous version of model Eq. 1, which is known to
support waves advecting a low number of _localised spikes_ , or having a non-
localised profile [31, 12, 14, 32, 65]. The travelling waves of interest to
us, however, have a localised profile, and advect a _large_ number of spikes,
such as the one presented in Fig. 1(a). These structures are not accessible
with the current techniques, hence we develop here analytical and numerical
tools to construct them. We define particular type of solutions, which retain
a fixed number of spikes in time; this class of solutions is sufficiently
general to incorporate travelling waves with an arbitrary, finite, number of
spikes, and small perturbations to them. We introduce the _voltage mapping_ ,
a new operator which formalises an idea previously used in the literature for
spiking [31, 12, 14, 32, 65, 3] and non-spiking networks [1, 35, 13, 21]. The
voltage mapping is based on level sets describing firing events, and it allows
efficient travelling wave constructions and stability computations.

Using the voltage mapping, we construct numerically waves with more than $200$
concurrent spikes. These waves are spatially localised, linearly unstable, and
coexist with the trivial (laminar) state (see Fig. 1(a)). As in the turbulence
analogy, the waves contain features of the bump attractor: they pack a
seemingly arbitrary number of spikes within the width of a bump attractor, and
they advect them at an arbitrarily slow speed, depending on $\beta$, and on
the number of carried spikes. As in the fluid-dynamical analogy, waves are
disconnected from the laminar state. Owing to the intrinsic non-smoothness of
the network, the waves emerge primarily at grazing points (as opposed to
saddle-node bifurcations seen in the fluid-dynamical analogy), which are
however seen in certain parameter regimes. In addition, we present numerical
evidence that the dynamics of the bump attractor displays echoes of the
unstable waves which, as in the fluid-dynamics analogy, form building blocks
for the localised structure. Also, the characteristic wandering of the bump
attractor, whose excursions become more prominent as $\beta$ increases, is
supported by this purely _deterministic_ system.

The paper is structured as follows: in Section 2 we introduce the discrete
model, characterise it as a non-smooth threshold network, and present
numerical simulations of bumps and waves; in Section 3 we introduce the
continuum model, the voltage mapping, and the construction of travelling
waves; in Section 4 we discuss travelling stability, we present numerical
results in Section 5, and we conclude in Section 6.

##  2 Coherent structures in the discrete model

We begin by introducing the discrete model by Laing and Chow [54]. We
characterise it as a piecewise-linear dynamical system, and we show numerical
simulations of coherent structures. An important difference from the work by
Laing and Chow is that we consider a _deterministic_ model, which we call the
Discrete Integrate-and-Fire Model (DIFM). We remark that the neurons
considered here, taken in isolation, are in an excitable regime, that is, they
exhibit an all-or-none response, based on the input they receive. This is
considerably different from the so-called oscillatory regime, in which
neurons, when decoupled from the network, display oscillations [82, 14, 62,
40].

###  2.1 Description of the DIFM

The DIFM is a spatially-extended system of $n$ identical _integrate-and-fire_
neurons, posed on $\mathbb{S}=\mathbb{R}/2L\mathbb{Z}$, that is, a ring of
period $2L$. Neurons are indexed using the set
$\mathbb{N}_{n}=\\{1,\ldots,n\\}$ and occupy the discrete, evenly spaced nodes
$x_{i}=-L+2iL/n\in\mathbb{S}$, for $i\in\mathbb{N}_{n}$. Neurons are coupled
via their synaptic connections, which are modelled by a continuous, bounded,
even and exponentially decaying function $w\colon\mathbb{S}\to\mathbb{R}$: the
strength of the connections from the $k$th to the $i$th neuron depends solely
on the distance $|x_{i}-x_{k}|$, measured around the ring, hence we write it
as $W_{ik}=w(x_{i}-x_{k})$, for all $i,k\in\mathbb{N}_{n}$.

To the $i$th neuron is associated a real-valued time-dependent voltage
function $v_{i}(t)$, and the coherent structures of interest are generated
when voltages $\\{v_{i}\\}$ attain a threshold value (when neurons _fire_).
The DIFM is formally written as follows:

(2) $\displaystyle\dot{v}_{i}(t)$
$\displaystyle=I_{i}(t)-v_{i}(t)+\frac{2L}{n}\sum_{k\in\mathbb{N}_{n}}\sum_{j\in\mathbb{N}}W_{ik}\alpha(t-\tau_{k}^{j})-\sum_{j\in\mathbb{N}}\delta(t-\tau^{j}_{i}),$
$\displaystyle i\in\mathbb{N}_{n},$ (3) $\displaystyle v_{i}(0)$
$\displaystyle=v_{0i},$ $\displaystyle i\in\mathbb{N}_{n}.$

At time $\tau_{i}^{j}$, when the voltage $v_{i}$ reaches the value $1$ from
below for the $j$th time, a firing event occurs; a more precise definition of
these _spiking times_ will be given below. The formal evolution equation Eq. 2
expresses the modelling assumption that, when a neuron fires, its voltage is
instantaneously reset to $0$ (hence the Dirac delta), and a so-called _post-
synaptic current_ is received by all other neurons in the network, with
intensity proportional to the strength of the synaptic connections. The time-
evolution of this current is modelled via the _post-synaptic function_
$\alpha(t)=p(t)H(t)$, expressed as the product of a continuous potential
function $p$ and the Heaviside function $H$, hence the post-synaptic current
is zero before a spike.

In this paper we present concrete calculations for

(4) $\alpha(t)=\beta\exp(-\beta t)H(t),\qquad
w(x)=a_{1}\exp(-b_{1}|x|)-a_{2}\exp(-b_{2}|x|),$

with $\beta,a_{1},a_{2},b_{1},b_{2}>0$, albeit the analytical and numerical
framework presented below is valid for more generic choices, subject to
general assumptions which will be made precise in Section 3.2. The function
$\alpha$ models exponentially-decaying currents with rate $-\beta$ and initial
value $\beta$, hence the limit $\beta\to\infty$ approximates instantaneous
currents. Currents with an exponential rise and decay are also used in
literature. The synaptic coupling function $w$ is chosen so that connections
are positive (_excitatory_) on the lengthscale $1/b_{1}$, and negative
(_inhibitory_) on the lengthscale $1/b_{2}$.

In addition to the post-synaptic current, neurons are subject to an external
stimulus $I_{i}(t)$. In certain time simulations, coherent structures will be
elicited with the application of a transient, heterogeneous stimulus of the
form

(5) $I_{i}(t)=I+d_{1}H(t-\tau_{\textrm{ext}})/\cosh(d_{2}x_{i}),\quad
i\in\mathbb{N}_{n}.$

Our investigation, however, concerns asymptotic states of the autonomous
homogeneous case $I_{i}(t)\equiv I$, hence one should assume $d_{1}=0$, unless
stated otherwise. A description of model parameters and their nominal values
can be found in Table 1.

Table 1: Parameter descriptions and nominal values. Parameter | Symbol | Value(s)
---|---|---
Number of neurons | $n$ | {80,500,1000,5000}
Domain half-width | $L$ | {1,3,4}
Synaptic efficacy and time scale | $\beta$ | [0,25]
Synaptic excitation coefficient | $a_{1}$ | 11
Synaptic inhibition coefficient | $a_{2}$ | 7
Synaptic excitation spatial scale | $b_{1}$ | 5
Synaptic inhibition spatial scale | $b_{2}$ | 3.5
Constant external input | $I$ | 0.9
Time-dependent external input duration | $\tau_{\textrm{ext}}$ | 2
Time-dependent external input strength | $d_{1}$ | {0,2}
Time-dependent external input spatial scale | $d_{2}$ | {10,12}

###  2.2 Event-driven DIFM

Laing and Chow studied and simulated a stochastic version of the DIFM, using
the Euler method and a first-order interpolation scheme to obtain the firing
times [54]. We use here a different approach: in preparation for our
analytical and numerical treatment of the problem, we write the formal model
Eqs. 2 to 3 as a system of $2n$ piecewise-linear ODEs. To this end we
introduce the synaptic input variables

(6)
$s_{i}(t)=\frac{2L}{n}\sum_{k\in\mathbb{N}_{n}}\sum_{j\in\mathbb{N}}W_{ik}\alpha(t-\tau_{k}^{j}),\qquad
i\in\mathbb{N}_{n}$

and combining Eq. 4 and Eq. 2 we obtain formally

$\begin{aligned}
\dot{v}_{i}(t)&=I_{i}(t)-v_{i}(t)+s_{i}(t)-\sum_{j\in\mathbb{N}}\delta(t-\tau^{j}_{i})\\\
\dot{s}_{i}(t)&=-\beta
s_{i}(t)+\frac{2L\beta}{n}\sum_{k\in\mathbb{N}_{n}}\sum_{j\in\mathbb{N}}W_{ik}\delta(t-\tau_{k}^{j})\end{aligned}\qquad
i\in\mathbb{N}_{n}.$

One way to define the associated non-smooth dynamical system is to express the
model as an impacting system, by partitioning the phase space
$\mathbb{R}^{2n}$ via a switching manifold, on which a reset map is prescribed
(see [27] and references therein for a discussion on non-smooth and impacting
systems). Here, we specify the dynamics so as to expose the firing times
$\\{\tau^{j}_{k}\\}$, as opposed to the switching manifold: this is natural in
the mathematical neuroscience context, and it prepares our analysis of the
continuum model. Since $\\{\tau^{j}_{k}\\}$ are the times at which orbits in
$\mathbb{R}^{2n}$ reach the switching manifold, a translation between the two
formalisms is possible.

Following these considerations, we set $\tau^{0}_{i}=0$ for all
$i\in\mathbb{N}_{n}$, introduce the notation $f({\mkern 2.0mu\cdot\mkern
2.0mu}^{\pm})=\lim_{\mu\to 0^{+}}f({\mkern 2.0mu\cdot\mkern 2.0mu}\pm\mu)$,
and define firing times as follows222 Note that $\\{\tau_{i}^{0}\\}_{i}$ are
not firing times, but auxiliary symbols for the definition of firing times Eq.
7. Indeed, since the sums in Eq. 2 run for $j\in\mathbb{N}$, the
$\\{\tau_{i}^{0}\\}_{i}$ are immaterial for the dynamics.

(7) $\tau_{i}^{j}=\inf\big{\\{}t\in\mathbb{R}\colon
t>\tau_{i}^{j-1},\;v_{i}(t^{-})=1,\;\dot{v}_{i}(t^{-})>0\big{\\}},\qquad
i\in\mathbb{N}_{n},\quad j\in\mathbb{N}.$

We arrange firing times in a monotonic increasing sequence
$\\{\tau_{i_{k}}^{j_{k}}\\}_{k=1}^{q}$ such that

(8)
$(0,T]=\bigcup_{k\in\mathbb{N}_{q+1}}\big{(}\tau^{j_{k-1}}_{i_{k-1}},\tau^{j_{k}}_{i_{k}}\big{]},\qquad
0=\tau^{j_{0}}_{i_{0}}<\tau^{j_{1}}_{i_{1}}\leq\ldots\leq\tau_{i_{q}}^{j_{q}}<\tau_{i_{q+1}}^{j_{q+1}}=T,$

for some time horizon $T>0$, and obtain the desired set of $2n$ piecewise-
linear ODEs

(9) $\dot{v}_{i}=I_{i}-v_{i}+s_{i},\quad\dot{s}_{i}=-\beta s_{i}\qquad
i\in\mathbb{N}_{n},\qquad
t\in\bigcup_{k\in\mathbb{N}_{q+1}}\big{(}\tau^{j_{k-1}}_{i_{k-1}},\tau^{j_{k}}_{i_{k}}\big{]},$

with initial and reset conditions

(10) $\displaystyle v_{i}(0)$ $\displaystyle=v_{0i},$ $\displaystyle s_{i}(0)$
$\displaystyle=s_{0i},$ $\displaystyle i\in\mathbb{N}_{n},$ (11)
$\displaystyle v_{i_{k}}(\tau^{j_{k}\,+}_{i_{k}})$ $\displaystyle=0,$
$\displaystyle s_{l}(\tau^{j_{k}\,+}_{i_{k}})$
$\displaystyle=s_{l}(\tau^{j_{k}\,-}_{i_{k}})+\frac{2L\beta}{n}W_{li_{k}},$
$\displaystyle l\in\mathbb{N}_{n},$ $\displaystyle k\in\mathbb{N}_{q},$

respectively. Henceforth, we refer to the non-smooth dynamical system Eqs. 7
to 11 with connectivity function $w$ given by Eq. 4 and stimulus Eq. 5 as the
_event-driven DIFM_ or simply _DIFM_ , that is, we view this model as a
substitute for the formal system Eqs. 2 to 3.

Even though the firing-time notation may seem cumbersome at first, the
evolution of the DIFM is remarkably simple: Equation 9 states that between two
consecutive firing times, neurons evolve independently, subject to a linear
ODE; a solution in closed form can be written in terms of exponential
functions, parametrised by the firing times. Constructing a solution amounts
to determining firing times (impacts with the switching manifold), as is
customary in piecewise-linear systems. This aspect will be a recurring theme
in the sections analysing travelling waves in the continuum model.

Figure 2: Bump attractors obtained via direct numerical simulation of the DIFM
Eqs. 7 to 11 with external input Eq. 5 and connectivity function $w$ as in Eq.
4. We visualise the network voltage (centre) and synaptic current (right) as
functions of space and time and, in the inset (left), a raster plot of the
firing events. Parameters as in Table 1 with $n=80$, $d_{1}=2$ $d_{2}=10$. The
network’s synaptic time scale is $\beta=1$ (a) and $\beta=3.5$ (b),
respectively. A localised coherent structure is visible in (a), which wanders
when $\beta$ is increased. We remark that the system under consideration is
deterministic.

In simulations of the DIFM, we time step Equation 9 rather than using its
analytic solution. We use an explicit adaptive 4-5th order Runge-Kutta pair
with continuous output, and detect events (compute firing times) by root-
finding [28, 76]. The simulation stops at each firing event and is restarted
after the reset conditions Eq. 11 are applied. Simulating the event-driven
DIFM instead of Eq. 2 allows us to compute firing times accurately, and to
evolve the system without storing in memory or truncating the synaptic input
sums in Eq. 2.

###  2.3 Coherent structures in the DIFM

Figure 3: Stable coexistent waves obtained via direct numerical simulation of
DIFM Eqs. 7 to 11 with external input Eq. 5 and connectivity function $w$ as
in Eq. 4. Parameters as in Table 1 with $n=80$, $\beta=4.5$ for both (a) and
(b), but different initial stimuli: (a) $d_{1}=0.4$, $d_{2}=12$, (b)
$d_{1}=2$, $d_{2}=10$. Depending on the transient stimulus the model displays:
(a) a wave propagating with positive speed, in which pairs of neurons fire
asynchronously, but at short times from each other; (b) a similar structure
involving a quartet of neurons. Coexisting structures with variable number of
firing neurons have also been found (not shown). The spatial profiles indicate
that neurons reach threshold (dashed red line) one at a time within a pair (a)
or two at a time within a quartet (b).

The DIFM supports standing and travelling localised structures, as in the
stochastic setting [54]. Bumps form robustly when we prescribe homogeneous
initial conditions333Typically we set $v_{0i}=u\in(0,1),s_{0i}=0$, for
$i\in\mathbb{N}_{n}$, but the coherent structures discussed in the paper can
also be found with random, independent and identically distributed initial
voltages, for instance $v_{0i}\sim\mathcal{U}([0,1])$, where $\mathcal{U}$ is
the uniform distribution. with a short transient stimulus (Equation 5 with
$\tau_{\textrm{ext}}=2$). Since $I_{i}(t)\equiv I$ for all
$t>\tau_{\textrm{ext}}$, the structures observed over long-time intervals are
solutions to a homogeneous, non-autonomous problem.

As seen in Fig. 2, the bump wanders when $\beta$ is increased. In passing, we
note that this phenomenon is not due to stochastic effects, as studied in
other contexts [50, 48, 3], because the DIFM is deterministic. For
sufficiently large $\beta$, the system exhibits stable travelling structures:
in Fig. 3 we show two coexisting waves, found for $\beta=4.5$ upon varying
slightly the width $d_{1}$ and intensity $d_{2}$ of the transient stimulus. In
each case we plot the voltage and synaptic profiles, and associated raster
plots. We notice different firing patterns in the waves, involving $2$ and $4$
firings, respectively: the wave with $2$ firings travels faster, and its
voltage and synaptic profiles are narrower. We found coexisting waves with a
greater number of firings and progressively lower speed, whose existence and
bifurcation structure will be at the core of the following sections.

###  2.4 Remarks about coherent structures in the DIFM

The patterns presented so far are found in the DIFM with a finite number of
neurons $(n=80)$. At first sight, the raster plots of the waves seem to
indicate that neurons fire simultaneously in pairs (Fig. 3(a)) or quartets
(Fig. 3(b)) as the structure travels across the network. A closer inspection
of the instantaneous profiles $v_{i}(t)$ reveal that this is not the case, as
the threshold (red dashed line) is attained by a single neuron in Fig. 3(a),
and by two neurons in Fig. 3(b): neurons in a raster pair fire alternatively
at short times from each other, whereas a quartet displays a more complex
firing pattern.

Hence, for finite $n$, the propagating structures displayed in Fig. 3 are not
strictly travelling waves, in the sense that the profile is not stationary in
the comoving frame; their dynamics is that of saltatory waves [22, 86, 3]. The
saltatory nature of the waves, however, is an effect of the network size: as
we increase $n$, the amplitude of temporal oscillations in the comoving frame
scales as $O(n^{-1})$, and the spatio-temporal profile converges to the one of
a travelling wave as $n\to\infty$.

In addition, the structure in Fig. 2(a) is not a bump, in the sense that it is
not a spatially heterogeneous steady state of the DIFM, because the pattern is
sustained by firing events. Indeed, the only equilibrium supported by the DIFM
is the homogeneous state $v_{i}(t)\equiv I$, $s_{i}(t)\equiv 0$,
$i\in\mathbb{N}_{n}$, which is linearly stable for all values of $\beta$, as
can be deduced by inspecting system Eq. 9.

It will be shown below that we can gain insight into the structure shown in
Fig. 2(a) (and its wandering) by constructing travelling waves and
investigating their stability in a continuum version of the DIFM.

##  3 Travelling waves in the continuum model

As stated in Section 2, the profiles $\\{v_{i}(t)\\}_{i}$ and
$\\{s_{i}(t)\\}_{i}$ in Fig. 3 behave like travelling wave solutions as
$n\to\infty$. Motivated by this observation, we study travelling waves in a
continuum, translation-invariant version of the DIFM: we set $d_{1}=0$ in the
stimulus Eq. 5, consider a continuum spatial domain, and pose the model on
$\mathbb{R}$ as opposed to $\mathbb{S}$, obtaining

(12)
$\begin{split}\partial_{t}v(x,t)=-v(x,t)+I&+\sum_{j\in\mathbb{N}}\int_{-\infty}^{\infty}w(x-y)\alpha\big{(}t-\tau_{j}(y)\big{)}\,dy\\\
&-\sum_{j\in\mathbb{N}}\delta\big{(}t-\tau_{j}(x)\big{)},\qquad(x,t)\in\mathbb{R}\times\mathbb{R}.\end{split}$

The formal evolution equation presented above, which we henceforth call the
_continuous integrate-and-fire model_ (CIFM), has been proposed and studied by
several authors in the mathematical neuroscience literature [31, 42, 10, 12,
65, 64]. In the CIFM, firing-time functions $\tau_{j}(x)$ indicate that the
neural patch at position $x$ fires for the $j$th time, and replace the
discrete model’s firing times $\tau_{k}^{j}$.444The index $j$ is used as a
superscript in the firing times, but for notational convenience we use it as a
subscript in the firing functions, so that
$\tau_{j}(x_{k})\approx\tau_{k}^{j}$. A graph of the firing functions replaces
the raster plot in the discrete model, so that a travelling wave in the CIFM
corresponding to the $n\to\infty$ limit of the structure in Fig. 3(a), for
instance, will involve $2$ linear firing functions $\tau_{1}$, $\tau_{2}$,
with $\tau_{1}(x)<\tau_{2}(x)$ for all $x\in\mathbb{R}$.

The existence of travelling waves solutions in Eq. 12 with a single spike has
been studied by Ermentrout [31] who presented various scalings of the
wavespeed as a function of control parameters. A general formalism for the
construction and linear stability analysis of _wavetrains_ (spatially-periodic
travelling solutions) was introduced and analysed by Bressloff [12], who
derived results in terms of Fourier series expansions. The construction of
travelling waves with multiple spikes was later studied by Oşan and coworkers
[64], albeit stability for these states was not presented and computations
were limited to a few spikes, for purely excitatory connectivity kernels. The
common thread in the past literature on this topic is the idea that travelling
wave construction and stability analysis rely entirely on knowledge of the
firing function $\tau_{j}$ (as in the DIFM, with firing times). A similar
approach has been used effectively in Wilson-Cowan-Amari neural field
equations, where it is often called _interfacial dynamics_ (see [Amari1977a]
for the first study of this type, [20] for a recent review, and [23, 38],
amongst others, for examples of spatio-temporal patterns analysis).

Here we present a new treatment of travelling wave solutions that draws from
this idea; we introduce an operator, that we call the _voltage mapping_ , with
the following aims: (i) Expressing a mapping between firing functions and
solution profiles, with the view of replacing the formal evolution equation
Eq. 12 for travelling waves with $m$ spikes (where $m$ is arbitrary). (ii)
Finding conditions for the linear stability of these waves. (iii) Using root-
finding algorithms to compute travelling waves and study their linear
stability. We will relate to existing literature in our discussion.

###  3.1 Notation

Before analysing solutions to the CIFM, we discuss the notation used in this
section. For fixed $m\in\mathbb{N}$, we use $|{\mkern 2.0mu\cdot\mkern
2.0mu}|_{\infty}$ to denote, with a little abuse of notation, the
$\infty$-norm on $\mathbb{C}^{m}$. We will use $|{\mkern 2.0mu\cdot\mkern
2.0mu}|$ for the standard modulus in $\mathbb{C}$. We denote by $C(X,Y)$ the
set of continuous functions from $X$ to $Y$, and use $C(X)$ when
$Y=\mathbb{R}$. We denote by $B(X)$ the set of real-valued bounded functions
defined on a set $X$, and by $BC(X)$ the set of real-valued, bounded and
continuous functions defined on $X$, respectively; both spaces are endowed
with the supremum norm $\|{\mkern 2.0mu\cdot\mkern 2.0mu}\|_{\infty}$. We
denote by $L^{1}(\mathbb{R})$ the pace of Lebesgue-integrable functions
defined on $\mathbb{R}$. Further, for a positive number $\eta$, we define the
exponentially weighted space

$L^{1}_{\eta}(\mathbb{R})=\Big{\\{}u\colon\mathbb{R}\to\mathbb{R}\colon\|u\|_{L^{1}_{\eta}}=\int_{\mathbb{R}}e^{\eta
x}|u(x)|\,dx<\infty\Big{\\}},$

which is a Banach space equipped with the norm $\|{\mkern 2.0mu\cdot\mkern
2.0mu}\|_{L^{1}_{\eta}}$, and the following space of exponentially growing
functions

$C_{\eta}(\mathbb{R},\mathbb{C}^{m})=\Big{\\{}u\in
C(\mathbb{R},\mathbb{C}^{m})\colon\|u\|_{C_{m,\eta}}=\sup_{x\in\mathbb{R}}e^{-\eta|x|}\,|u(x)|_{\infty}<\infty\Big{\\}},$

which is a Banach space equipped with the norm $\|{\mkern 2.0mu\cdot\mkern
2.0mu}\|_{C_{m,\eta}}$.

###  3.2 Characterisation of solutions to the CIFM via the voltage mapping

We begin by discussing in what sense a voltage function $v$ satisfies the CIFM
formal evolution equation Eq. 12. While we eschew the definition of the CIFM
as a dynamical system on a Banach space (a characterisation that is currently
unavailable in the literature), we note that progress can be made for voltage
profiles with a constant and finite number of spikes for $t\in\mathbb{R}$.
This class of solutions is sufficiently large to treat travelling waves, and
small perturbations to them.

We make a few assumptions on the network coupling, and we restrict the type of
firing functions and solutions of interest, as follows: [Coupling functions]
The connectivity kernel $w$ is an even function in $C(\mathbb{R})\cap
L^{1}_{\eta}(\mathbb{R})$, for some $\eta>0$. The post-synaptic function
$\alpha\colon\mathbb{R}\to\mathbb{R}_{\geq 0}$ can be written as
$\alpha(t)=p(t)H(t)$, where $H$ is the Heaviside function, and
$p\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$ is a bounded and everywhere
differentiable Lipschitz function, hence $p,p^{\prime}\in B(\mathbb{R})$.

######  Definition 3.1 ($m$-spike CIFM solution).

Let $m\in\mathbb{N}$ and $I\in\mathbb{R}$. A function
$v_{m}\colon\mathbb{R}^{2}\to\mathbb{R}$ is an $m$-spike CIFM solution if
there exists $\tau=(\tau_{1},\ldots,\tau_{m})\in C(\mathbb{R},\mathbb{R}^{m})$
such that $\tau_{1}<\ldots<\tau_{m}$ on $\mathbb{R}$ and

(13) $\displaystyle\begin{aligned}
v_{m}(x,t)=I&+\sum_{j\in\mathbb{N}_{m}}\int_{-\infty}^{t}\int_{-\infty}^{\infty}\\!\\!\\!\exp(z-t)w(x-y)\alpha(z-\tau_{j}(y))\,dy\,dz\\\
&-\sum_{j\in\mathbb{N}_{m}}\exp(\tau_{j}(x)-t)H(t-\tau_{j}(x)),\qquad(x,t)\in\mathbb{R}^{2}\end{aligned}$
(14) $\displaystyle v_{m}(x,t)=1,\qquad(x,t)\in\mathbb{F}_{\tau},$ (15)
$\displaystyle
v_{m}(x,t)<1,\qquad(x,t)\in\mathbb{R}^{2}\setminus\mathbb{F}_{\tau},$

where

$\mathbb{F}_{\tau}=\bigcup_{j\in\mathbb{N}_{m}}\\{(x,t)\in\mathbb{R}^{2}\colon
t=\tau_{j}(x)\\}.$

We call $\tau$ and $\mathbb{F}_{\tau}$ the firing functions and the firing set
of $v_{m}$, respectively.

The definition above specifies how we interpret solutions to Eq. 12, and is
composed of three ingredients: (i) Equation 13, which derives from integrating
Eq. 12 on $(-\infty,t)$, and expresses a mapping between the set of $m$ firing
functions $\tau$ and the voltage profile; (ii) System Eq. 14, which couples
the firing functions by imposing the threshold crossings; (iii) A further
condition on $v_{m}$, ensuring that the solution has exactly $m$ spikes,
attained at the firing set; this is necessary because, as we shall see below,
it is possible to find a set of $m$ functions $\tau$ satisfying Equations 13
to 14, but exhibiting a number of threshold crossings greater than $m$.

We now aim to characterise $m$-spike CIFM solutions by means of a voltage
mapping, which can be conveniently linearised around a firing set, and is a
key tool to construct waves and analyse their stability. Inspecting Eq. 13 we
note that the voltage profile features two contributions, one from the
(synaptic) coupling functions $w$ and $\alpha$, and one from reset conditions.
This observation leads to the following definitions:

######  Definition 3.2 (Synaptic and reset operators).

Let $u:\mathbb{R}\to\mathbb{R}$. We define the synaptic operator, $S$, and the
reset operator, $R$, by

(16) $\displaystyle(Su)(x,t)$
$\displaystyle=\int_{-\infty}^{t}\int_{-\infty}^{\infty}\exp(z-t)w(x-y)\alpha(z-u(y))\,dy\,dz,$
$\displaystyle(x,t)\in\mathbb{R}^{2},$ (17) $\displaystyle(Ru)(x,t)$
$\displaystyle=-\exp(u(x)-t)H(t-u(x)),$ $\displaystyle(x,t)\in\mathbb{R}^{2}.$

These operators map univariate functions, such as a firing function, to
bivariate functions, contributing to the spatio-temporal voltage profile. The
following lemma shows that the synaptic contribution is a continuous function
on the plane, hence discontinuities in the voltage come through the reset
operator, as expected.

######  Lemma 3.3.

If Section 3.2 holds, then for the operators $S$, $R$ in Definition 3.2 we
have $S\colon C(\mathbb{R})\to BC(\mathbb{R}^{2})$ and $R\colon
C(\mathbb{R})\to B(\mathbb{R}^{2})$, respectively.

######  Proof 3.4.

The statement for $R$ is immediate, whereas proving the continuity of $Su$ on
$\mathbb{R}^{2}$ requires some estimates for the improper integral. A proof is
given in Appendix A.

We can now define the voltage mapping as follows, combining $S$ and $R$.

######  Definition 3.5 (Voltage mapping).

Let $m\in\mathbb{N}$, $I\in\mathbb{R}$ and $\tau\in
C(\mathbb{R},\mathbb{R}^{m})$. The $m$-spike voltage mapping, $V_{m}$, is the
operator defined as

(18) $V_{m}\tau=I+\sum_{j\in\mathbb{N}_{m}}(S\tau_{j}+R\tau_{j}),$

where $S$ and $R$ are given in Definition 3.2.

By construction, the voltage operator characterises $m$-spike CIFM solutions,
as the following proposition shows.

######  Proposition 3.6.

Let $m\in\mathbb{N}$, $I\in\mathbb{R}$. An $m$-spike CIFM solution exists if,
and only if, there exists $\tau\in C(\mathbb{R},\mathbb{R}^{m})$ such that

(19) $\displaystyle V_{m}\tau=1,$ $\displaystyle\text{in
$\mathbb{F}_{\tau}$},$ (20) $\displaystyle V_{m}\tau<1,$ in
$\mathbb{R}^{2}\setminus\mathbb{F}_{\tau}$

######  Proof 3.7.

The statement follows by setting $v_{m}(x,t)=(V_{m}\tau)(x,t)$ and applying
the definition of the voltage mapping, Equation 18.

In some cases it is useful to replace the threshold conditions Eqs. 14 and 19
by equivalent conditions involving left limits of the voltage function and
mapping, respectively, as specified by the following result.

######  Corollary 3.8 (Discontinuities of $v_{m}$).

Under the hypotheses of Lemma 3.3, $v_{m}=1$ in $\mathbb{F}_{\tau}$ if, and
only if, $\lim_{\mu\to
0^{+}}v_{m}(x,\tau_{i}(x)-\mu)=v_{m}(x,\tau_{i}(x)^{-})=1$ for all
$(i,x)\in\mathbb{N}_{m}\times\mathbb{R}$.

######  Proof 3.9.

Using Lemma 3.3 one can show that
$v_{m}(x,\tau_{i}(x))=v_{m}(x,\tau_{i}(x)^{-})$, as we derive in Appendix B.

Proposition 3.6 implies that the voltage of an $m$-spike solution can be
computed for any $(x,t)\in\mathbb{R}^{2}$ once the firing functions $\tau$ are
known. The spatio-temporal profile of an $m$-spike solution is determined
entirely by its firing functions. This aspect, which underlies the formal
evolution equation Eq. 12 and the literature which analyses it, is a key part
of what follows and, as we shall see below, it also suggests a natural way to
compute travelling waves, and determine their linear stability. A first step
in this direction is the definition of travelling waves via the voltage
mapping.

###  3.3 Travelling waves with m-spikes ($\textrm{TW}_{m}$)

Following Proposition 3.6, we can capture travelling waves with $m$ spikes
($\textrm{TW}_{m}$) using the voltage mapping, and a set of parallel firing
functions. Henceforth, we will assume without loss of generality that the
propagating speed of the wave is positive: for any wave with $c>0$, there
exists a wave with speed $-c$, and the wave profiles related by the
transformation $x\to-x$.

######  Definition 3.10 ($\textrm{TW}_{m}$).

Let $m\in\mathbb{N}$, $c>0$, and let $T\in\mathbb{R}^{m}$ with
$T_{1}<\cdots<T_{m}$. A travelling wave with $m$ spikes ($\textrm{TW}_{m}$),
speed $c$, and coarse variables $(c,T)$ is an $m$-spike CIFM solution with
firing functions $\\{\tau_{j}(x)=x/c+T_{j}\\}_{j\in\mathbb{N}_{m}}$.

To each travelling wave solution is associated a travelling wave profile which
is advected with propagation speed $c$. From Proposition 3.6 we expect this
profile to be determined entirely by the firing functions, as confirmed in the
following result.

######  Proposition 3.11 ($\textrm{TW}_{m}$ profile).

A $\textrm{TW}_{m}$ with speed $c$ satisfies
$(V_{m}\tau)(x,t)=\nu_{m}(ct-x;c,T)$, and its $(c,T)$-dependent travelling
wave profile $\nu_{m}$ is given by

(21) $\displaystyle\nu_{m}(\xi;c,T)=I$
$\displaystyle-\sum_{j\in\mathbb{N}_{m}}\exp\bigg{(}-\frac{\xi-
cT_{j}}{c}\bigg{)}H\bigg{(}\frac{\xi-cT_{j}}{c}\bigg{)}$
$\displaystyle+\frac{1}{c}\sum_{j\in\mathbb{N}_{m}}\int_{-\infty}^{\xi}\exp\bigg{(}\frac{z-\xi}{c}\bigg{)}\int_{0}^{\infty}w(y-z+cT_{j})p(y/c)\,dy\,dz.$

######  Proof 3.12.

The statement follows from evaluating the operators $S$ and $R$ at the firing
functions $\\{\tau_{j}(x)=x/c+T_{j}\\}_{j\in\mathbb{N}_{m}}$, and operating a
change of variables. See Appendix C.

Proposition 3.11 shows that the travelling wave profile is completely
determined by the vector $(c,T)\in\mathbb{R}_{>0}\times\mathbb{R}^{m}$, that
is, $(c,T)$ is a vector of coarse variables for the travelling wave. In the
discrete model we introduced an auxiliary spatially-extended variable for the
model, the synaptic input $\\{s_{i}(t)\\}_{i}$ defined in Eq. 6. In the
continuum model, the corresponding variable is the function
$s_{m}(x,t)=\sum_{j\in\mathbb{N}_{m}}(S\tau_{j})(x,t)$, which in a
$\textrm{TW}_{m}$ satisfies $s_{m}(x,t)=\sigma_{m}(ct-x;c,T)$, with

(22)
$\sigma_{m}(\xi;c,T)=\frac{1}{c}\sum_{j=1}^{m}\int_{0}^{\infty}w(y-\xi+cT_{j})p(y/c)\,dy.$

###  3.4 Travelling wave construction

Figure 4: Wave profiles for (a) $\textrm{TW}_{5}$ and (b) $\textrm{TW}_{20}$
obtained by solving Problem 3.13 for $m=5$ and $m=20$, respectively, and then
subsitituting $(c,T_{1},\dots,T_{m})$ into the expression for voltage profile
Eq. 21 and synaptic profile Eq. 22. The profile $\nu$ is computable at any
$\xi\in\mathbb{R}$, here we plot it using an arbitrary grid in the intervals
(a) $[-0.5,1]$ and (b) $[-1.5,2]$. Parameters as in Table 1 with (a)
$\beta=4.5$ and (b) $\beta=7.7$.

Proposition 3.11 suggests a simple way to compute a $\textrm{TW}_{m}$, by
determining its $m+1$ _coarse variables_ $(c,T)$, as a solution to the
following _coarse problem_ :

######  Problem 3.13 (Computation of $\textrm{TW}_{m}$).

Find $(c,T)\in\mathbb{R}_{>0}\times\mathbb{R}^{m}$ such that
$T_{1}<\cdots<T_{m}$ and

(23) $\displaystyle T_{1}$ $\displaystyle=0,$ (24)
$\displaystyle\nu_{m}(cT_{i}^{-};c,T)$ $\displaystyle=1,\qquad\text{for
$i\in\mathbb{N}_{m}$},$ (25) $\displaystyle\nu_{m}(\xi;c,T)$
$\displaystyle<1,\qquad\text{on
$\mathbb{R}\setminus\cup_{j\in\mathbb{N}_{m}}\\{cT_{j}^{-}\\}$}.$

Equation 24 of the coarse problem imposes that the travelling wave profile
crosses the threshold $1$ when $\xi\to cT^{-}_{j}$. As expected, if $\nu_{m}$
is a travelling wave profile, then so is $\nu_{m}(\xi+\xi_{0})$ for any
$\xi_{0}\in\mathbb{R}$; Equation 23 fixes the phase of the travelling wave, by
imposing that the profile crosses threshold as $\xi\to 0^{-}$.

If $m=1$, Equations 23 to 24 of the coarse problem reduce to a compatibility
condition for the speed $c$,

$c\int_{-\infty}^{0}\int_{0}^{\infty}\exp(s)w\big{(}c(y-s)\big{)}p(y)\,dy\,ds=I-1,$

which implicitly defines an existence curve for $\textrm{TW}_{1}$ in the
($c$,$I$)-plane. This result is in agreement with what was found in [64, 31].
Existence curves in other parameters are also possible, and are at the core of
the numerical bifurcation analysis presented in detail in the sections below.

For $m>1$, the coarse problem must be solved numerically. A simple solution
strategy is to find a candidate solution using Newton’s method for the system
of $m+1$ transcendental equations Eqs. 23 to 24, with $\nu_{m}$ given by
Proposition 3.11, and with initial guesses estimated from direct simulation of
the discrete model with large $n$, or from a previously computed coarse
vector. The candidate solution can then be evaluated at arbitrary
$\xi\in\mathbb{R}$, hence it is accepted if Eq. 24 holds on a spatial grid
covering $[-L,L]\subset\mathbb{R}$, with $L\gg 1$. In passing, we note that
this procedure is considerably cheaper than a standard travelling wave
computation for PDEs, which requires the solution of a boundary value problem,
and hence a discretisation of differential operators on $\mathbb{R}$.
Depending on the particular choice of $\alpha$ and $w$, the profile $\nu_{m}$
is either written in closed form, as is the case for the choices Eq. 4, or
approximated using standard quadrature rules.

A concrete calculation is presented in Fig. 4, where we show travelling wave
profiles and speeds of a $\textrm{TW}_{5}$ and a $\textrm{TW}_{20}$. In
passing, we note that the synaptic profile of a $\textrm{TW}_{m}$ at a given
time is similar to a bump, but displays modulations at the core (visible in
Fig. 4), as predicted by the Heaviside switches in Eq. 22. Travelling waves
with a large number of spikes, such as these ones, have not been accessible to
date.

######  Remark 3.14.

Fig. 4 shows that profiles with $\nu_{m}(cT_{j}^{-})=1$ propagate with
positive speed, and this does not contradict the numerical simulations in Fig.
3, where solutions profiles with $v_{m}(x,\tau_{j}(x)^{-})=1$ propagate with
negative speed. This is a consequence of choosing $\xi=ct-x$ (as in [64]),
hence initial conditions for the time simulations are obtained by reflecting
$\nu_{m}$ about the $y$ axis, since $v_{m}(x,0)=\nu_{m}(-x)$.

##  4 Wave Stability

Figure 5: (a)-(b): Examples illustrating the destabilisation of a
$\textrm{TW}_{3}$ solution. A time simulation of the DIFM is initialised using
wave profiles obtained solving Problem 3.13 for $m=3$ at (a) $\beta=17$ and
(b) $\beta=17.5$. Parameters as in Table 1, domain half-width $L=4$ and
network size $n=1000$. The firing functions $\\{\tau_{j}\\}$ are plotted for
reference. Oscillatory perturbations to the firing functions do not decrease
with time, hence the wave is unstable. The dynamics leads to stable (a)
$\textrm{TW}_{2}$ and (b) $\textrm{TW}_{1}$ solutions. (c): Perturbations
$\tau+\varphi$ to the firing functions $\tau$ of a $\textrm{TW}_{m}$. At $t=0$
each firing function $\tau_{i}$ is perturbed by an amount
$\varphi_{i}(-cT_{i})$. A $\textrm{TW}_{m}$ is linearly stable if
$\varphi_{i}(-cT_{i})$ being small implies that $\varphi_{i}(x)$ stays small
for all $x\in(-cT_{i},\infty)$ and $i\in\mathbb{N}_{m}$ (see Definition 4.4).

The time simulations in Section 2 demonstrate that, for sufficiently large
values of $\beta$, travelling waves with a variable number of spikes coexist
and are stable. It is natural to ask whether these waves destabilise as
$\beta$, or any other control parameter of the model, is varied. An example of
a prototypical wave instability is presented in Fig. 5 for $\textrm{TW}_{3}$:
a travelling wave is computed solving Problem 3.13, and this solution is used
as initial condition for a DIFM simulation with $n=1000$ neurons. For
sufficiently large $\beta$, the wave in unstable, as exemplified by the raster
plots in Fig. 5(a)–(b), in that the firing function never return to the ones
of a $\textrm{TW}_{3}$.

Inspecting Fig. 5 we observe that the firing set of the solution is composed
of 3 disjoint curves. Initially, these curves are close to the ones of a
$\textrm{TW}_{3}$, from which they depart progressively. Ultimately, some
firing functions terminate, and the dynamics displays an attracting
$\textrm{TW}_{2}$ or $\textrm{TW}_{1}$. Capturing the transitions from a
$\textrm{TW}_{m}$ to a travelling wave with fewer spikes is a nontrivial task.
Firstly, classifying these transitions implies studying the _nonlinear
stability_ of the wave. Secondly, the characterisation we have provided of the
CIFM solutions, Definition 3.1, assumes that the number of spikes is constant
at all times: accounting for changes in numbers of the firing functions would
require a redefinition of the model, starting from the formal evolution
equation Eq. 12, and patching together solutions segments with variable
numbers of spikes.

The voltage mapping, however, opens up the possibility of studying the _linear
stability_ of $\textrm{TW}_{m}$: the spatio-temporal voltage profile of an
$m$-spike solution is determined by its firing functions, $\tau$, via Eq. 18;
small perturbations $\tau+\varepsilon\varphi$ to $\tau$, induce small
perturbations to the spatio-temporal profile, and we expect that a suitable
linearisation of the voltage mapping carries information concerning the
asymptotic behaviour of these perturbations.

The main aim of this section is to formalise the concept of linear stability
for the problem under consideration, and to provide an algorithm for
$\textrm{TW}_{m}$ linear stability computations. We begin by showing that if
two distinct $m$-spike solutions have firing functions $\tau$ and
$\tau+\varphi$, respectively, then the perturbations $\varphi$ satisfy a
linear equation to leading order. The following lemma also specifies
admissible perturbations, namely $\varphi$ are in the Banach space
$C_{\eta}(\mathbb{R},\mathbb{R}^{m})$: perturbations are allowed to grow
exponentially as $|x|\to\infty$, at a rate at most equal to $\eta$, which
bounds the decay rate of the connectivity kernel function $w\in
L^{1}_{\eta}(\mathbb{R})$.

######  Lemma 4.1 (Linearised voltage mapping operator).

Assume Section 3.2, and let $(c,T)$ be the coarse variables of a
$\textrm{TW}_{m}$ with firing functions $\tau$. Further, let $L$ be the linear
operator defined by
$L\varphi=\big{(}(L\varphi)_{1},\ldots,(L\varphi)_{m}\big{)}$, where

$(L\varphi)_{i}=\sum_{j\in\mathbb{N}_{m}}(\varphi_{i}-\varphi_{j})1_{j<i}+\int_{cT_{ji}}^{\infty}e^{-y/c}w(y)\psi_{ij}(y)\big{[}\varphi_{i}-\varphi_{j}({\mkern
2.0mu\cdot\mkern 2.0mu}-y)\big{]}\,dy,\quad i\in\mathbb{N}_{m},$

with coefficients $T_{ij}$ and functions $\psi_{ij}$ given by

$T_{ij}=T_{i}-T_{j},\qquad\psi_{ij}\colon[cT_{ij},\infty)\to\mathbb{R},\quad
y\mapsto p(0)+\int_{0}^{y/c-T_{ji}}e^{s}p^{\prime}(s)\,ds,\qquad
i,j\in\mathbb{N}_{n},$

respectively. The following statements hold:

  1. 1.

$L$ is a bounded operator from $C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ to
itself.

  2. 2.

Let $0<\varepsilon\ll 1$ and $\varphi\in C_{\eta}(\mathbb{R},\mathbb{R}^{m})$.
If $\tau+\varepsilon\varphi$ are firing functions of an $m$-spike CIFM
solution (a perturbation of the $\textrm{TW}_{m}$), then

(26) $0=L\varphi+O(\varepsilon)\qquad\text{in $\mathbb{R}$}$

######  Proof 4.2.

To prove that $L\colon C_{\eta}(\mathbb{R},\mathbb{C}^{m})\to
C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ is well defined and bounded, one shows
that for any $\varphi\in C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ there exists a
positive constant $\kappa_{m,\eta}$ such that
$\|L\varphi\|_{C_{m,\eta}}\leq\kappa_{m,\eta}\|\varphi\|_{C_{m,\eta}}$. A
proof is obtained by obtaining preliminary estimates on $\psi_{ij}$, and
$|(L\varphi)(x)|_{\infty}$, and relies on the decaying properties of $w$,
stated in Section 3.2. Part 2 is proved by setting
$u=\tau+\varepsilon\varphi\in C_{\eta}(\mathbb{R},\mathbb{R}^{m})$,
linearising the mapping $u\mapsto V_{m}u-1|_{\mathbb{F}_{u}}$ around $\tau$,
and exploiting properties of $m$-spikes solutions outlined in the Section 3.2.
A proof is given in Appendix D.

######  Remark 4.3.

Note that the operator $L$ depends on the coarse variables $(c,T)$, albeit we
omit this dependence for notational simplicity.

We are now ready to define linear stability for a $\textrm{TW}_{m}$, which we
adapt from [12]. Intuitively, we compare the firing set $\mathbb{F}_{\tau}$ of
a $\textrm{TW}_{m}$ with the firing set $\mathbb{F}_{\tau+\varphi}$ of a
perturbed $m$-spike solution with $\|\varphi\|_{C_{m,\eta}}\ll 1$, for which
$\varphi$ satisfy Eq. 26 to leading order. If the sets $\mathbb{F}_{\tau}$ and
$\mathbb{F}_{\tau+\varphi}$ are close around $t=0$ and remain close for all
positive times, we deem the wave linearly stable. With reference to Fig. 5(c),
we observe that, when $\textrm{TW}_{m}$ crosses the axis $t=0$, each one of
its firing functions $\tau_{i}$ is perturbed by an amount
$\varphi_{i}(-cT_{i})$. Roughly speaking, a $\textrm{TW}_{m}$ is linearly
stable if $\varphi_{i}(-cT_{i})$ being small implies that $\varphi_{i}(x)$
stays small for all $x\in(-cT_{i},\infty)$ and $i\in\mathbb{N}_{m}$. If a wave
is linearly stable and all $\varphi_{i}$ decay to $0$ as $x\to\infty$ we say
that the wave is asymptotically linearly stable. More precisely:

######  Definition 4.4 (Linear stability of $\textrm{TW}_{m}$).

A $\textrm{TW}_{m}$ with coarse varaibles $(c,T)$ is linearly stable to
perturbations $\varphi$ if $\varphi\in\ker L$, and for each $\varepsilon>0$
there exists $\delta=\delta(\varepsilon)>0$, such that if
$|\varphi_{i}(-cT_{i})|<\delta$, then $|\varphi_{i}(x)|<\varepsilon$ for all
$(i,x)\in\mathbb{N}_{m}\times(-cT_{i},\infty)$.

A $\textrm{TW}_{m}$ is asymptotically linearly stable to perturbations
$\varphi$ if is linearly stable to perturbations $\varphi$ and
$|\varphi(x)|_{\infty}\to 0$ as $x\to\infty$.

We have seen that a $\textrm{TW}_{m}$ can be constructed by solving a
nonlinear problem in the unknowns $(c,T)$. The following lemma, which is the
central result of this section, establishes that linear stability of a
$\textrm{TW}_{m}$ with respect to exponential perturbations of the firing
functions can also be determined by finding roots of a $(c,T)$-dependent,
complex-valued function.

######  Lemma 4.5 ($\textrm{TW}_{m}$ stability).

Assume Section 3.2, let $(c,T)$ be coarse variables of a $\textrm{TW}_{m}$,
and let $\mathbb{D}_{a,b}=\\{z\in\mathbb{C}\colon a\leq\operatorname{Re}z\leq
b\\}$. Further, let $E$ be the complex-valued function

(27) $E\colon\mathbb{D}_{-\eta,\eta}\to\mathbb{C},\quad z\mapsto\det[D-M(z)],$

where $M\in\mathbb{C}^{m\times m}$,
$D=\operatorname{diag}(D_{1},\ldots,D_{m})\in\mathbb{R}^{m\times m}$, are the
matrices with elements

$M_{ij}(z)=e^{T_{ji}}\bigg{[}1_{j<i}+\int_{cT_{ji}}^{\infty}e^{-(z+1/c)y}w(y)\psi_{ij}(y)\,dy\bigg{]},\qquad
D_{i}=\sum_{k\in\mathbb{N}_{m}}M_{ik}(0),$

respectively, then:

  1. 1.

If $\lambda$ is a root of $E$, then its complex conjugate $\lambda^{*}$ is
also a root of $E$, and there exists a nonzero $\Phi\in\ker[D-M(\lambda)]$
such that $\Phi e^{\lambda x},\Phi^{*}e^{\lambda^{*}x}\in\ker{L}$, where $L$
is defined as in Lemma 4.1.

  2. 2.

$E$ has a root at $0$. $\textrm{TW}_{m}$ is linearly stable (but not
asymptotically linearly stable) to perturbations $\varphi\colon x\mapsto\kappa
v$, where $\kappa\in\mathbb{R}\setminus\\{0\\}$ and
$v=(1,\ldots,1)\in\mathbb{R}^{m}$.

  3. 3.

If $\lambda$ is a root of $E$ in
$\mathbb{D}_{-\eta,0}\setminus\mathrm{i}\mkern 1.0mu\mathbb{R}$, then
$\textrm{TW}_{m}$ is linearly asymptotically stable to perturbations $\Phi
e^{\lambda x}+\Phi^{*}e^{\lambda^{*}x}$.

######  Proof 4.6.

A proof of the statement is given in Appendix E.

Lemma 4.5 provides a link between exponential perturbations to the firing
times of a $\textrm{TW}_{m}$ and zeroes of the function $E$ in the strip
$\mathbb{D}_{-\eta,\eta}\subset\mathbb{C}$. The function $E$ depends on
$(c,T)$ via the entries of the matrices $D,M$, and can be evaluated
numerically at each point $z\in\mathbb{D}_{-\eta,\eta}$.

In PDEs, linear stability of a travelling wave is determined by the spectrum
of a linear operator, which contains a $0$ eigenvalue corresponding to a
translational perturbation mode. Part 2 of Lemma 4.5 provides an analogous
result for a $\textrm{TW}_{m}$, which is linearly stable, but not
asymptotically linearly stable (therefore neutrally stable), to perturbations
that shift the firing functions homogeneously. Part 3 of Lemma 4.5 suggests
that a $\textrm{TW}_{m}$ is stable if all nonzero roots of $E$ have strictly
negative real parts. Initial guesses for the roots can be obtained by plotting
$0$-level sets of the function $E$, for fixed $(c,T)$.

##  5 Bifurcation structure of travelling waves

Figure 6: (a) Branch of $\textrm{TW}_{3}$ solutions in the parameter $\beta$,
using $c$ as solution measure. The branch originates at a grazing point $G$,
illustrated by the profile in (b). As $\beta$ increases, three pairs of
complex conjugate roots of $E$ (see Equation 27) cross the imaginary axis at
the oscillatory (Hopf) bifurcation points $\textrm{HB}_{1}$,
$\textrm{HB}_{2}$, $\textrm{HB}_{3}$. Panels (c) and (d) show selected roots
of $E$, before and after $\textrm{HB}_{1}$, at $\beta=10$ and $16$,
respectively. (e)–(g) Raster plots for time simulations of the DIFM with
$n=500$ and domain half-width $L=3$, initialised from solutions to Problem
3.13 at $\beta=2.17$, $10$, and $16$, respectively. The simulations show the
dynamics of the model for $\beta<\beta_{\textrm{G}}$ (where a
$\textrm{TW}_{3}$ does not exist in the continuum limit), for
$\beta\in(\beta_{\textrm{G}},\beta_{\textrm{HB}_{1}})$ (where
$\textrm{TW}_{3}$ is stable according to the analysis in (c)), and for
$\beta>\beta_{\textrm{HB}_{1}}$ (where $\textrm{TW}_{3}$ is unstable to
oscillatory perturbations, as predicted in (d)). Parameters as in Table 1,
with $d_{1}=0$.

The pseudo-arclength continuation routines developed in [67, 2] have been used
to compute solutions to Problem 3.13, continue waves in parameter space, and
investigate their stability. A $\textrm{TW}_{m}$ is constructed by solving
Problem 3.13 in the corase variables
$(c,T)\in\mathbb{R}_{>0}\times\mathbb{R}^{m}$, which is sufficient to
reconstruct the wave profile Eq. 21, and the corresponding synaptic profile
Eq. 22; in addition, starting from a solution to Problem 3.13, the linear
asymptotic stability of a $\textrm{TW}_{m}$ is determined by finding roots of
the $(c,T)$-dependent nonlinear function $E$ defined in Eq. 27.

Fig. 6 shows the bifurcation structure of $\textrm{TW}_{3}$, which is common
to most travelling waves found in the model. The simulations in Section 2
suggest to take the synaptic timescale parameter $\beta$ as the principal
continuation parameter. We use the wavespeed $c$ as solution measure. A branch
of solutions originates from a grazing point (G, see below for a more detailed
explanation) and it is initially stable, before destabilising at a sequence of
oscillatory bifurcations ($\textrm{HB}_{1}$–$\textrm{HB}_{3}$), as seen in
Fig. 6(a). In passing, we note that there exists a second, fully unstable,
branch of $\textrm{TW}_{3}$ solutions characterised by a slower speed and a
smaller width. This branch, which we omit from the bifurcation diagrams for
simplicity, also originates at a grazing point.

###  5.1 Grazing points

In a wide region of parameter space, branches of $\textrm{TW}_{m}$ solutions
originate at a grazing point $\beta=\beta_{\textrm{G}}$, as seen in Fig.
6(a)–(b) for $\textrm{TW}_{3}$555Note that $\beta_{\textrm{G}}$ depend on $m$,
but we omit this dependence to simplify notation. The same is true for other
quantities in the paper such as $c$ and $T_{\textrm{G}}$, for instance.. At a
grazing point the $\textrm{TW}_{m}$ profile crosses threshold $m$ times, and
attains the threshold tangentially at a further spatial location,
$cT_{\textrm{G}}$, as shown in Fig. 6(b). This tangency exists at the critical
value $\beta=\beta_{\textrm{G}}$, signalling a non-smooth transition and a
branch termination. For $\beta>\beta_{\textrm{G}}$ we observe profiles with
exactly $m$ threshold crossings (a branch of $\textrm{TW}_{m}$ solutions).
These profiles exhibit a further local maximum, which is strictly less than
$1$ by construction, at a point $\xi_{\textrm{max}}>cT_{m}$. As
$\beta\to\beta^{+}_{\textrm{G}}$ we observe $\xi_{\textrm{max}}\to
cT_{\textrm{G}}^{+}$ and $\nu(\xi_{\textrm{max}})\to 1^{-}$, until the
threshold is reached at $\beta=\beta_{\textrm{G}}$, where the tangency
originates.

For $\beta<\beta_{\textrm{G}}$ we find solutions to the nonlinear problem Eqs.
23 to 24 for which $V_{m}\tau<1$ in a bounded interval of $\mathbb{R}$. Since
these states violate the condition Eq. 25, they do not correspond to
$\textrm{TW}_{m}$ solutions, and we disregard them (the branch terminates at
$\beta_{G}$). We note, however, that in a neighbourhood of
$\beta_{\textrm{G}}$ there exist branches of travelling wave solutions with
different number of threshold crossings (as it will be shown below).

We found grazing points for every $\textrm{TW}_{m}$ with $2\leq m\leq 230$,
for the parameters in Table 1 with $d_{1}=0$. We observe that for
$\beta<\beta_{\textrm{G}}$ the system evolves towards a DIFM bump attractor
(see Fig. 6(e)). Understanding the origin of this transition is the subject of
the following sections.

Grazing points are found generically as a secondary control parameter is
varied, and $2$-parameter continuations of grazing points can be obtained
numerically, by freeing one parameter and imposing tangency of the wave
profile at one additional point (see Problem F.1 in Appendix F).

###  5.2 Oscillatory bifurcations

Along the $\textrm{TW}_{m}$ branch, we compute and monitor the roots of $E$
with the largest real part. Fig. 6(c)-(d) show examples for $\textrm{TW}_{3}$
at $\beta=10$ and $\beta=16$ respectively. At $\beta=10$ we observe a root at
$0$, as expected, and other roots with small negative real part: the wave is
therefore linearly asymptotically stable to firing-threshold perturbations
$x\mapsto\Phi e^{\lambda x}+\Phi^{*}e^{\lambda^{*}x}$, with $E(\lambda)=0$ and
$\Phi\in\ker[D-M(\lambda)]$, as confirmed via simulation in Fig. 6(f). In
contrast, there exists a pair of unstable complex conjugate roots for the
solution at $\beta=16$, indicating an oscillatory (Hopf) instability, which is
also confirmed by direct simulation, in Fig. 6(g): after the initial
oscillatory instability, the system destabilises to a $\textrm{TW}_{2}$. It
should be noted that, in other regions of parameter space and for simulations
with different network sizes, we observed a $\textrm{TW}_{3}$ destabilise to a
$\textrm{TW}_{1}$ or the homogeneous steady state.

We expect that branches of periodically modulated $\textrm{TW}_{m}$ solutions
(which are also supported by neural fields [33, 24]) emerge from each of the
Hopf bifurcations reported in Fig. 6(a). We note that we could not find stable
structures of this type via direct simulations near onset and, while it is
possible to extend our framework to continue such periodic states, we did not
pursue this strategy here.

As shown in Fig. 6(a), the $\textrm{TW}_{3}$ branch undergoes a sequence of
Hopf bifurcations $\\{\textrm{HB}_{i}\\}_{i}$: our stability analysis shows
several pairs of complex conjugate roots progressively crossing the imaginary
axis as $\beta$ increases: the computation in Fig. 6(d), for instance, is for
a solution at $\beta\in(\beta_{\textrm{HB}_{1}},\beta_{\textrm{HB}_{2}})$. We
have verified numerically (not shown) that the firing functions of spatio-
temporal DIFM solutions in this region of parameter behave as predicted by the
leading eigenvalues in Fig. 6(d), that is, they feature two dominant
oscillatory modes: one stable, and one unstable. Similarly to grazing points,
Hopf bifurcations can be continued in a secondary parameter (see Problem F.2
in Appendix F).

###  5.3 Nested branches of travelling waves

Figure 7: Bifurcation structure of $\textrm{TW}_{m}$ branches for
$m=1,\ldots,160$ in the parameter $\beta$. (a): For $m\geq 3$, branches are
similar to the one shown in Fig. 6(a). As $m$ increases, the waves become
slower and their stability region narrower. The shaded area in (a) is enlarged
in (b): the inset shows selected branches for $m=2,\ldots,160$; oscillatory
instabilities occur within the red segments, and the branches with $m\geq 57$
are fully unstable (solid grey lines). We used here the same data as in Fig.
1, but we present it in terms of $c$, not $\Delta$. Parameters as in Table 1,
with $d_{1}=0$.

We computed branches of $\textrm{TW}_{m}$ solutions for increasing values of
$m$, as reported in Fig. 7(a), using DIFM simulations as initial guesses. In
Fig. 1 waves were represented by their width, whereas here we use the
propagation speed $c$. In the region of parameter space explored in the DIFM
model, branches with $m\geq 2$ feature a grazing point for low $\beta$, and
branches with $m\geq 3$ display sequences of Hopf Bifurcations, following the
scenario already discussed in Fig. 6(a). In this region, the $\textrm{TW}_{1}$
branch has a distinct behaviour, featuring a saddle node point in place of a
grazing point. For each $\textrm{TW}_{m}$ branch terminating at a grazing
point, there is a corresponding slow unstable branch originating at a
different grazing point: in Fig. 7(a) this behaviour is exemplified by
plotting the fully unstable slow $\textrm{TW}_{5}$ branch (the branch with
slowest waves in the figure), but is omitted for all other branches. The two
$\textrm{TW}_{5}$ branches should be understood as a “broken saddle-node”. The
bifurcation structure of Fig. 7(a), valid for the CIFM, supports numerical
simulations of the DIFM, in which a $\textrm{TW}_{m}$ destabilises at
$\textrm{HB}_{1}$, and gives rise to a new travelling wave state,
$\textrm{TW}_{m^{\prime}}$ with $m^{\prime}<m$ (see for instance Figs. 5 and
6).

These _coexisting $\textrm{TW}_{m}$ branches_ are nested in a characteristic
fashion, so far unreported in the literature; the higher $m$, the slower the
wave, and the narrower the stable interval between $G$ and $\textrm{HB}_{1}$.
This structure is noteworthy: firstly, it is known that the speed of
$\textrm{TW}_{1}$ typically changes _as a secondary parameter is varied_ [31,
11, 12]; however, in networks with purely excitatory kernels, waves with
multiple threshold crossings coexist, and their speed does not depend strongly
on $m$ [31], which has been a principle reason for studying approximately and
analytically the only tractable case, $m=1$ [13, Section 5.4] (this scenario
is also confirmed by our calculations, see Fig. 11); secondly, it is known
that Hopf instabilities with purely excitatory connectivity kernel are
possible only if delays are present in the network [12].

The results in Fig. 7 have been obtained using a methodology that works for
arbitrary $m$, and on generic connectivity kernels. They show that, when
inhibition is present: (i) coexisting nested branches of $\textrm{TW}_{m}$
exist; (ii) the speed of such waves depends strongly on $m$, and in particular
it is possible to construct waves with arbitrarily small speed, by increasing
the number of spikes; (iii) oscillatory instabilities are present in models
without delays, for sufficiently large $m$ and/or sufficiently large $\beta$.
As we shall see, the latter aspect plays a role in understanding the so called
_bump attractor_.

Figure 8: (a) The quantities $c$ and $T_{m}$, evaluated at the grazing points
$\beta=\beta_{G}$, are $O(m^{-1})$ and $O(m)$, respectively. Since $T_{1}=0$
for all waves, the quantity $cT_{m}$ measures the wave width, and we expect
the sequence $\\{cT_{m}\\}_{m\in\mathbb{N}}$, the sequence of wave widths, to
converge to a fixed value as $m\to\infty$. (b) The solid grey line is the gain
function, proposed in [54] for a stationary asynchronous bump, red dots mark
the instantaneous firing rate for $\textrm{TW}_{230}$ at the grazing point,
computed according to the formula $(T_{i+1}-T_{i})^{-1}$ at position
$x=cT_{i}$.

###  5.4 The bump attractor

From the grazing point of $\textrm{TW}_{m}$, one can compute the grazing point
of $\textrm{TW}_{m+1}$. For instance, from the $\textrm{TW}_{3}$ grazing
profile in Fig. 6(b), we obtain $(c,T_{1},T_{2},T_{3},T_{G})$. A grazing point
can then be computed solving Problem F.1, and its solution can be used to
produce an initial guess $(c,T_{1},T_{2},T_{3},(T_{3}+T_{G})/2,T_{G})$ for a
grazing point of $\textrm{TW}_{4}$. Exploiting this iterative strategy, we
computed grazing points and branches for large values of $m$, obtaining the
diagram in Fig. 7(b), corresponding to the shaded area in Fig. 7(a).

The branches accumulate as $m$ increase, and for $m\geq 57$, they are fully
unstable for this parameter set. The diagrams provide evidence that there
exists unstable waves with arbitrarily many spikes and vanishingly small
speed. It seems therefore natural to postulate a relationship between these
waves and the structures found by Laing and Chow [54] (see also Figs. 2 and
6(e)).

We initially explore this relationship by inspecting the travelling wave
profiles: since each $\textrm{TW}_{m}$ branch with $m>1$ terminates at a
grazing point, we examine a few features of the profile at these points. The
leftmost spike of each wave occurs at $\xi_{1}=0$ by construction (see Problem
3.13), while its righmost spike is at $\xi_{m}=cT_{m}$, which is therefore a
proxy for the wave’s width666Recall that $c$ is also a function of $m$, but we
omit this dependence for ease of notation. Fig. 8(a) shows $c$ and $T_{m}$,
computed at the grazing points, as functions of $m$: we find $c=O(m)$,
$T_{m}=O(m^{-1})$, therefore we expect the sequence
$\\{\xi_{m}\\}_{m\in\mathbb{N}}$ to converge to a finite value $\xi_{*}$ as
$m\to\infty$. These data show that, as the wavespeed tends to zero, the
growing number of spikes are distributed in a fixed interval $[0,\xi_{*}]$
resembling therefore a stationary bump of width $\xi_{*}$. This conclusion is
also supported by an inspection of the distribution of spikes
$\\{\xi_{i,m}\\}_{i}$ as $m\to\infty$. Instead of looking directly at this
distribution, we plot the inverse inter-spike time $1/(T_{i+1}-T_{i})$, as a
function of $\xi_{i}$, for the grazing point of $\textrm{TW}_{230}$. Laing and
Chow [54] call this quantity the _gain function_ , and propose it as an
emerging firing rate function for the asynchronous bumps of the DIFM [54].
Importantly, Laing and Chow show that the gain function for asynchronous bumps
resembles the one of a neural network rate equation, but it is an emergent
property of the DIFM. Fig. 8(b) shows that gain function by Laing and Chow is
in excellent agreement with the one obtained for $\textrm{TW}_{230}$ (the
slowest wave we computed in the CIFM), confirming that, from a macroscopic
viewpoint, asynchronous bumps could be understood in the limit $m\to\infty$ of
$\textrm{TW}_{m}$.

Figure 9: (a): Mean instantaneous speed ($\bar{c}$ in Eq. 28, purple dots) and
interval estimators ($[\bar{c}-\sigma_{c},\bar{c}+\sigma_{c}]$ and
$[c_{\textrm{min}},c_{\textrm{max}}]$, dark and light purple shades,
respectively) in direct simulations of the DIFM, superimposed on
$\textrm{TW}_{m}$ branches of the CIFM (an inset of Fig. 7(b), which has been
reflected by the $c=0$ axis to signpost waves with negative speed). The bump
attractor is characterised by $\bar{c}\approx 0$, and speed fluctuations which
grow with $\beta$, while echoing the oscillatory instabilities of the branch
(red segments). (b): exemplary solutions in (a) displaying an initial
advection, followed by a bump attractor (1,2) or a stable wave (3). (c):
Histograms of the instantaneous speed in selected time intervals, indicated by
the lateral coloured bars in (b). We observe transitions through weakly
unstable waves (orange interval in (3), blue intervals in (1, 2)),
corresponding to sharp peaks in the histograms.

We have further investigated the bump attractor state, in relation to the
$\textrm{TW}_{m}$: the analysis of the continuum model, in the region of
parameter space where the bump attractor is observed, predicts the coexistence
of the trivial attracting solution $v(x,t)\equiv I$, with arbitrarily slow,
unstable waves, whose spatial profile approximates the one of a bump. We
simulated the DIFM with $n=5000$, initialising the model from an unstable
travelling wave of the CIFM, $\textrm{TW}_{105}$, and estimated the
instantaneous speed $c(t)$ of the numerical DIFM solution at $q$ time points
$\\{t_{k}\colon k\in\mathbb{N}_{q}\\}$, using a level set of the synaptic
profile and finite differences, as follows:

$z(t)=\max\\{x\in\mathbb{S}:s(x,t)=0.1\\},\qquad
c_{k}=(z(t_{k})-z(t_{k-1}))/(t_{k}-t_{k-1}),\qquad k\in\mathbb{N}_{q}.$

A CIFM travelling wave solution corresponds to a constant $c$: when the DCLM
solution displays a wave for large $n$, the sequence $\\{c_{k}\\}_{k}$
converges to a constant value, if one disregards small oscillations due to the
finite $n$, and which vanish as $n\to\infty$. On the other hand we expect that
no differentiable function $c(t)$ exists for a bump attractor, albeit some
information may be contained in the mean, $\overline{c}$, standard deviation,
$\sigma_{c}$, and extrema, $c_{\textrm{min}}$, $c_{\textrm{max}}$ of the
_deterministic_ scalar $c_{k}$

(28)
$\bar{c}=\frac{1}{q}\sum_{k\in\mathbb{N}_{q}}c_{k},\quad\sigma^{2}_{c}=\frac{1}{q-1}\sum_{k\in\mathbb{N}_{q}}(c_{k}-\bar{c})^{2},\quad
c_{\textrm{min}}=\min_{k\in\mathbb{N}_{q}}c_{k},\quad
c_{\textrm{max}}=\max_{k\in\mathbb{N}_{q}}c_{k},$

which we computed for various values of $\beta$, and superimposed on the
bifurcation diagram of the CIFM model, in Fig. 9(a): we plot $\bar{c}$ (purple
dots) and two interval estimators, $[\bar{c}-\sigma_{c},\bar{c}+\sigma_{c}]$
(dark purple shade) and $[c_{\textrm{min}},c_{\textrm{max}}]$ (light purple
shade); we recall that the CFIM admits branches of waves with positive and
negative speed, both plotted in the figure. The bump attractor is
characterised by a zero-mean speed, and interval estimators which grow with
$\beta$, indicating that the instantaneous speed undergoes larger fluctuations
as the bump meanders. We observe that the interval estimators echo the onset
of oscillations (red segments) in the $\textrm{TW}_{m}$ branches, as $\beta$
increases.

This behaviour is robust, for low and medium values of $\beta$, albeit the
fine details of the dynamics depend on initial conditions. Figure 9(b) shows 3
examples whose estimated speeds appear also in Fig. 9(a). The space-time plots
display an initial advection, followed by a bump attractor, or a stable
travelling wave. To highlight the transitions, we computed histograms of
$c_{k}$ in selected time intervals, indicated by blue, orange, yellow, and
purple bars in Fig. 9(b). The histograms provide numerical evidence of
transitions near weakly unstable waves (orange intervals/histograms in example
3, blue intervals/histograms in examples 1 and 2), which manifest themselves
via sharp peaks.

###  5.5 Composite waves

In addition to the waves studied thus far, we found by direct simulation waves
whose firing functions are split into well-separated groups, that is, firing
functions in the same group are closer to each another than they are to those
in other groups, see Fig. 10. We call these structures composite waves, as
they may be formed via the interaction of travelling waves with various
numbers of spikes. As in other non-smooth dynamical systems [44], we expect
that these solutions have discontinuities that are rearranged with respect to
a $\textrm{TW}_{m}$.

For illustrative purposes, we denote a composite waves with $k\in\mathbb{N}$
groups by $\textrm{TW}_{m_{1}}$ \+ …+ $\textrm{TW}_{m_{k}}$, where
$\\{m_{i}\\}_{i=1}^{k}$ is a sequence of positive integers specifying the
number of spikes in each group. There are constraints for the groups, dictated
by dynamical considerations: for instance a $\textrm{TW}_{1}$ \+
$\textrm{TW}_{3}$ cannot exist, because a $\textrm{TW}_{1}$, taken in
isolation, is faster than a $\textrm{TW}_{3}$. The construction of asymptotic
profiles and computation of linear stability for composite waves follow in the
same way as defined in Sect. 3 and Sect. 4.

In Fig. 10(a), we show a selection of of composite waves near the
$\textrm{TW}_{3}$ branch. Roughly speaking, the wave profile along each
depicted branch comprises a $\textrm{TW}_{3}$ as its leading group, followed
by two additional spike group that collectively form a compound satisfying the
travelling wave conditions (e.g. branch 1 combines a $\textrm{TW}_{3}$, a
$\textrm{TW}_{2}$ and a $\textrm{TW}_{1}$). The branches of composite waves
are separate from each other and from the previously computed
$\textrm{TW}_{m}$ branches in Fig. 7, however, all branches possess a
bifurcation structure similar to the one of the $\textrm{TW}_{m}$ discussed in
the past section. Moreover, we see that the magnitude of the speed of the
composite wave is bounded above by the magnitude of the speed of the group at
the leading edge of the wave (the slowest wave, $\textrm{TW}_{3}$ in this
case).

Figure 10: (a) Bifurcation diagram of selected composite waves. The red curve
is a $\textrm{TW}_{3}$ branch, as computed in Fig. 6. The blue curves are
branches of composite waves, featuring an approximate $\textrm{TW}_{3}$ at the
front of the wave. The composite waves are slightly slower than
$\textrm{TW}_{3}$. The diagram shows selected profiles at the first
oscillatory bifurcation points. (b) Examples of composite waves obtained via
collisions of multi-spike waves. (c) Collisions between $m$-spike propagating
structures and wandering bumps generate composite waves (left) or bump
repulsion (right), depending on initial conditions. Simulations in panels
(b)–(c) have a lattice spacing of $\Delta x=2L/n=0.01$.

Direct numerical simulation highlights that composite waves can be formed from
the interaction of multi-spike waves as shown in the left panel of Fig. 10(b).
Here we choose an initial condition with well separated $\textrm{TW}_{1}$,
$\textrm{TW}_{2}$ and $\textrm{TW}_{3}$ profiles. Initially, these separated
structures travel with different speeds ($\textrm{TW}_{1}$ being the fastest
and $\textrm{TW}_{3}$ the slowest, in line with what was found in Fig. 7(a)).
After a transient, the waves come closer and form a compound (the composite
wave), with a common intermediate speed. The dynamics of composite waves
depend greatly on the initial conditions: in the right panel of Fig. 10(b), we
see that an initial condition in which a $\textrm{TW}_{1}$ lies between
another $\textrm{TW}_{1}$ and a $\textrm{TW}_{3}$ leads to the extinction of
the intermediate wave resulting in a composite wave with a total of $4$
spikes.

Composite waves also result from the collision between waves and wandering
bumps (Fig. 10(c), left panel). Here, we see a transition of two bump states
into a composite wave that is compounded with a pre-existing
$\textrm{TW}_{5}$. The interaction with the $\textrm{TW}_{5}$ causes the left-
most bump to visit the branches of travelling wave solutions whereupon the
combined state settles on a stable $\textrm{TW}_{5}$ \+ $\textrm{TW}_{9}$
branch. This process is repeated for the right-most bump, giving rise to an
overall $\textrm{TW}_{5}$ \+ $\textrm{TW}_{9}$ \+ $\textrm{TW}_{12}$ wave. In
the right panel of the Fig. 10(c), we see that the same kind of collision can
instead result in the wave packet transitioning to a wandering bump itself,
highlighting the dependence of the formation of composite waves on initial
conditions. In this scenario, the bump state does not visit a stable
travelling wave branch and so only transiently adopts a weakly unstable wave
profile before returning to a bump attractor state.

##  6 Conclusions

We have provided evidence that the relationship between bump attractors and
travelling waves in a classical network of excitable, leaky integrate-and-fire
neurons bears strong similarities to the one between complex spatiotemporal
patterns and waves at the onset of pipe turbulence. We made analytical and
numerical progress in the construction and stability analysis of travelling
waves with a large number of localised spikes, and gained access to their
intricate bifurcation structure. This step was essential, because such waves
advect, at low speed, localised patterns that resemble the bump attractor
core. It should be noted that the waves we computed are only a subset of the
ones supported by the model.

As we completed the present paper, a recent publication [55] reported the
existence of waves with vanishingly small speed, and discontinuous profiles,
in networks of theta neurons, which can be cast as spiking networks with a
polynomial ODE of quadratic type. A natural question arises as to whether the
fluid-dynamical analogy applies in that and other network models. The level-
set approach used in the present paper was particularly effective because one
can define $m$-spike waves starting from mild solutions to the formal
evolution equation Eq. 12, and derive a relatively simple expression for the
wave profile Eq. 21. While this approach may be harder to carry out in more
detailed spiking models, the general idea of a relationship between localised
waves and bumps in spiking networks could be investigated, by direct
simulations, in more realistic networks (spiking or not).

An important open question concerns the definition of Eq. 12 and, more
generally, of spatially-continuous spiking networks, as dynamical systems
posed on function spaces. This problem has been circumvented here by defining
a suitable class of solutions, introducing the voltage mapping, and then
providing proofs of its relevance to the construction and stability of
multiple-spike waves. We believe that a full dynamical-systems
characterisation of similar models will be a key ingredient to uncover further
links between localised waves and bumps in complex, spatially-extended
threshold networks.

## Acknowledgments

We are grateful to Stephen Coombes, Predrag Cvitanović, Gregory Faye, Joel
Feinstein, John Gibson, Joost Hulshof, and Edgar Knobloch. for insightful
discussions.

##  Appendix A Proof of Lemma 3.3

######  Lemma A.1 (Re-statement of Lemma 3.3).

If Section 3.2 holds, then for the operators $S$, $R$ in Definition 3.2 we
have $S\colon C(\mathbb{R})\to BC(\mathbb{R}^{2})$ and $R\colon
C(\mathbb{R})\to B(\mathbb{R}^{2})$, respectively.

######  Proof A.2.

Fix $u\in C(\mathbb{R})$. The real-valued function $z\mapsto\exp(-z)H(z)$ is
bounded in $\mathbb{R}$, hence $Ru\in B(\mathbb{R}^{2})$. To prove the result
on $S$, we define the functions

$\psi(y,t)=\int_{\infty}^{t}e^{z-t}\alpha(z-u(y))\,dz,\qquad
s(x,t)=(Su)(x,t)=\int_{-\infty}^{\infty}w(x-y)\psi(y,t)\,dy,$

and set $K_{\alpha}=\|\alpha\|_{\infty}$, $K_{w}=\|w\|_{L^{1}(\mathbb{R})}$,
whose existence is guaranteed by Section 3.2. The function $s$ is bounded on
$\mathbb{R}^{2}$ because $|\psi|\leq K_{\alpha}$ on $\mathbb{R}^{2}$, hence

$|s(x,t)|\leq K_{\alpha}\int_{-\infty}^{\infty}|w(x-y)|\,dy\leq
K_{\alpha}K_{w}.$

In order to prove the continuity of $s$, it is useful to first show that
$\psi$ is continuous in $t$ on $\mathbb{R}$, uniformly in $y$. This claim is
proved by noting that, for any $y,t,\tau\in\mathbb{R}$ with $t<\tau$ we have

$\displaystyle|\psi(y,\tau)-\psi(y,t)|$
$\displaystyle=\bigg{|}\int_{-\infty}^{t}[e^{z-\tau}-e^{z-t}]\alpha(z-u(y))\,dz+\int_{t}^{\tau}e^{z-\tau}\alpha(z-u(y))\,dz\bigg{|}$
$\displaystyle\leq
K_{\alpha}\int_{-\infty}^{t}e^{z-t}-e^{z-\tau}\,dz+K_{\alpha}(\tau-t)$
$\displaystyle=K_{\alpha}\big{(}1-e^{-(\tau-t)}+\tau-t\big{)},$

which, combined with a similar argument for $\tau<t$, leads to

(29) $|\psi(y,\tau)-\psi(y,t)|\leq
K_{\alpha}\big{(}1-e^{-|\tau-t|}+|\tau-t|\big{)},\qquad\text{for all
$y,t,\tau\in\mathbb{R}$}.$

We prove the continuity of $s$ by showing that $|s(\xi,\tau)-s(x,t)|\to 0$ as
$(\xi,\tau)\to(x,t)$. We consider the following inequality

(30) $|s(\xi,\tau)-s(x,t)|\leq|s(\xi,\tau)-s(x,\tau)|+|s(x,\tau)-s(x,t)|,$

and we note that the second term in the right-hand side of Eq. 30 can be made
arbitrarily small as $(\xi,\tau)\to(x,t)$, owing to Eq. 29. Therefore it
suffices to show that the first term in the right-hand side of Eq. 30 can also
be made arbitrarily small as $(\xi,\tau)\to(x,t)$, that is, we must show that
for any $\varepsilon>0$ there exists $\delta>0$ such that
$|s(\xi,\tau)-s(x,\tau)|<\varepsilon$. To prove this statement, we use the
boundedness of $\psi$ and a change of variables in the integral to obtain the
estimate

$\displaystyle|s(\xi,\tau)-s(x,\tau)|$
$\displaystyle\leq\int_{-\infty}^{\infty}|w(y-\xi)-w(y-x)|\,|\psi(y,\tau)|\,dy$
$\displaystyle\leq K_{\alpha}\int_{-\infty}^{\infty}|w(y+x-\xi)-w(y)|\,dy,$

therefore, for any $X>0$ we have

$\displaystyle|s(\xi,\tau)-s(x,\tau)|$ $\displaystyle\leq
K_{\alpha}\int_{-\infty}^{-X}|w(y+x-\xi)-w(y)|\,dy$
$\displaystyle+K_{\alpha}\int_{-X}^{X}|w(y+x-\xi)-w(y)|\,dy$
$\displaystyle+K_{\alpha}\int_{X}^{\infty}|w(y+x-\xi)-w(y)|\,dy:=K_{\alpha}(I_{1}+I_{2}+I_{3}).$

We bound $I_{3}$ as follows

$\displaystyle I_{3}$
$\displaystyle\leq\int_{X}^{\infty}|w(y+x-\xi)|\,dy+\int_{X}^{\infty}|w(y)|\,dy$
$\displaystyle=\int_{X+x-\xi}^{\infty}|w(y)|\,dy+\int_{X}^{\infty}|w(y)|\,dy\leq
2\int_{X-|x-\xi|}^{\infty}|w(y)|\,dy$

and a similar reasoning gives an identical bound for $I_{1}$,

$I_{1}\leq 2\int_{X-|x-\xi|}^{\infty}|w(y)|\,dy.$

We conclude that, for any $X>0$

$|s(\xi,\tau)-s(x,\tau)|\leq
4K_{\alpha}\int_{X-|x-\xi|}^{\infty}|w(y)|\,dy+K_{\alpha}\int_{-X}^{X}|w(y+x-\xi)-w(y)|\,dy.$

We now fix $\varepsilon,\delta_{1}>0$. Since $w\in L^{1}(\mathbb{R})$, we can
pick $X$ so that

$|x-\xi|<\delta_{1}\qquad\Rightarrow\qquad
4K_{\alpha}\int_{X-|x-\xi|}^{\infty}|w(y)|\,dy<\frac{\varepsilon}{2}.$

Furthermore, by continuity of $w$ there exists $\delta_{2}>0$ such that

$|x-\xi|<\delta_{2}\qquad\Rightarrow\qquad
K_{\alpha}\int_{-X}^{X}|w(y+x-\xi)-w(y)|\,dy\leq\frac{\varepsilon}{2},$

hence for any $\varepsilon>0$ there exists
$\delta=\min(\delta_{1},\delta_{2})>0$ such that
$|s(\xi,\tau)-s(x,\tau)|<\varepsilon$, which implies the continuity of $s$. We
conclude that $S\colon C(\mathbb{R})\to BC(\mathbb{R}^{2})$.

##  Appendix B Proof of Corollary 3.8

######  Corollary B.1 (Restatement of Corollary 3.8).

Under the hypotheses of Lemma 3.3, $v_{m}=1$ in $\mathbb{F}_{\tau}$ if, and
only if, $\lim_{\mu\to
0^{+}}v_{m}(x,\tau_{i}(x)-\mu)=v_{m}(x,\tau_{i}(x)^{-})=1$ for all
$(i,x)\in\mathbb{N}_{m}\times\mathbb{R}$.

######  Proof B.2.

The condition $v_{m}(x,t)=1$ for $(x,t)\in\mathbb{F}_{\tau}$ is equivalent to

$v_{m}(x,\tau_{i}(x))=1\qquad(i,x)\in\mathbb{N}_{m}\times\mathbb{R}.$

From the defintions of $S$ and $R$, and the continuity $S\tau_{i}$ on
$\mathbb{R}^{2}$ (see Lemma 3.3) we have, for all
$(i,x)\in\mathbb{N}_{m}\times\mathbb{R}$

$\displaystyle v_{m}(x,\tau_{i}(x))$
$\displaystyle=I+\sum_{j\in\mathbb{N}_{m}}(S\tau_{j})(x,\tau_{i}(x))-\sum_{j<i}\exp(\tau_{j}(x)-\tau_{i}(x))$
$\displaystyle=I+\sum_{j\in\mathbb{N}_{m}}\lim_{\mu\to 0^{+}}
(S\tau_{j})(x,\tau_{i}(x)-\mu)-\sum_{j<i}\exp(\tau_{j}(x)-\tau_{i}(x))$
$\displaystyle=I+\lim_{\mu\to
0^{+}}\sum_{j\in\mathbb{N}_{m}}\bigg{[}(S\tau_{j})(x,\tau_{i}(x)-\mu)-(R\tau_{j})(x,\tau_{i}(x)-\mu)\bigg{]}=v_{m}(x,\tau_{i}(x)^{-}),$

hence $v_{m}=1$ in $\mathbb{F}_{\tau}$ if, and only if,
$v_{m}(x,\tau_{i}(x)^{-})=1$ for all $(i,x)\in\mathbb{N}_{m}\times\mathbb{R}$.

##  Appendix C Proof of Proposition 3.11

######  Proposition C.1 ($\textrm{TW}_{m}$ profile).

A $\textrm{TW}_{m}$ with speed $c$ satisfies
$(V_{m}\tau)(x,t)=\nu_{m}(ct-x;c,T)$, and its $(c,T)$-dependent travelling
wave profile $\nu_{m}$ is given by

(31) $\displaystyle\nu_{m}(\xi;c,T)=I$
$\displaystyle-\sum_{j\in\mathbb{N}_{m}}\exp\bigg{(}-\frac{\xi-
cT_{j}}{c}\bigg{)}H\bigg{(}\frac{\xi-cT_{j}}{c}\bigg{)}$
$\displaystyle+\frac{1}{c}\sum_{j\in\mathbb{N}_{m}}\int_{-\infty}^{\xi}\exp\bigg{(}\frac{z-\xi}{c}\bigg{)}\int_{0}^{\infty}w(y-z+cT_{j})p(y/c)\,dy\,dz.$

######  Proof C.2.

We set $\tau_{j}(x)=x/c+T_{j}$ for $j\in\mathbb{N}_{n}$. The first sum in Eq.
31 is immediate, as

$(R\tau_{j})(x,t)=\exp\bigg{(}\frac{cT_{j}-ct+x}{c}\bigg{)}H\bigg{(}\frac{ct-
x-cT_{j}}{c}\bigg{)}$

The second sum is obtained as follows:

$\displaystyle(S\tau_{j})(x,t)$
$\displaystyle=\int_{-\infty}^{t}\int_{-\infty}^{\infty}\exp(s-t)w(x-y)\alpha(s-\tau_{j}(y))\,dy\,ds$
$\displaystyle=\int_{-\infty}^{t}\int_{0}^{\infty}\exp(s-t)w(x+y-cs+cT_{j})p(y/c)\,dy\,ds$
$\displaystyle=\frac{1}{c}\int_{-\infty}^{ct-x}\int_{0}^{\infty}\exp\bigg{(}\frac{s+x-ct}{c}\bigg{)}w(y-s+cT_{j})p(y/c)\,dy\,ds.$

##  Appendix D Proof of Lemma 4.1

######  Lemma D.1 (Restatement of Lemma 4.1).

Assume Section 3.2, and let $(c,T)$ be the coarse variables of a
$\textrm{TW}_{m}$ with firing functions $\tau$. Further, let $L$ be the linear
operator defined by
$L\varphi=\big{(}(L\varphi)_{1},\ldots,(L\varphi)_{m}\big{)}$, where

$(L\varphi)_{i}=\sum_{j\in\mathbb{N}_{m}}(\varphi_{i}-\varphi_{j})1_{j<i}+\int_{cT_{ji}}^{\infty}e^{-y/c}w(y)\psi_{ij}(y)\big{[}\varphi_{i}-\varphi_{j}({\mkern
2.0mu\cdot\mkern 2.0mu}-y)\big{]}\,dy,\quad i\in\mathbb{N}_{m},$

with coefficients $T_{ij}$ and functions $\psi_{ij}$ given by

$T_{ij}=T_{i}-T_{j},\qquad\psi_{ij}\colon[cT_{ij},\infty)\to\mathbb{R},\quad
y\mapsto p(0)+\int_{0}^{y/c-T_{ji}}e^{s}p^{\prime}(s)\,ds,\qquad
i,j\in\mathbb{N}_{n},$

respectively. The following statements hold:

  1. 1.

$L$ is a bounded operator from $C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ to
itself.

  2. 2.

Let $0<\varepsilon\ll 1$ and $\varphi\in C_{\eta}(\mathbb{R},\mathbb{R}^{m})$.
If $\tau+\varepsilon\varphi$ are firing functions of an $m$-spike CIFM
solution (a perturbation of the $\textrm{TW}_{m}$), then

(32) $0=L\varphi+O(\varepsilon)\qquad\text{in $\mathbb{R}$}$

######  Proof D.2.

_Part 1_. If $\varphi\in C_{\eta}(\mathbb{R},\mathbb{C}^{m})$, then
$\varphi\in C(\mathbb{R},\mathbb{C}^{m})$ and $L\varphi\in
C(\mathbb{R},\mathbb{C}^{m})$. We show that for any $\varphi\in
C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ there exists a positive constant
$\kappa_{m,\eta}$ such that
$\|L\varphi\|_{C_{m,\eta}}\leq\kappa_{m,\eta}\|\varphi\|_{C_{m,\eta}}$, which
implies that $L$ is a bounded operator from
$C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ to itself. We begin by estimating
$\psi_{ij}$: by Section 3.2 there exist constants $K_{p}$, $K_{p^{\prime}}$
such that

$\displaystyle|\psi_{ij}(y)|$
$\displaystyle\leq|p(0)|+\int_{0}^{y/c-T_{ji}}e^{z}|p^{\prime}(z)|\,dz$
$\displaystyle\leq
K_{p}+K_{p^{\prime}}\int_{0}^{y/c-T_{ji}}e^{z}\,dz=K_{p}+K_{p^{\prime}}(e^{y/c-T_{ji}}-1),$

therefore, introducing the constant
$K=2\max(K_{p},K_{p^{\prime}})\max_{i,j}e^{-T_{ji}}$,

(33) $e^{-y/c}|\psi_{ij}(y)|\leq
e^{-y/c}(K_{p}+K_{p^{\prime}}e^{y/c-T_{ji}})\leq K$

uniformly in
$(i,j,y)\in\mathbb{N}_{m}\times\mathbb{N}_{m}\times[cT_{ji},\infty)$.

We now fix $\varphi\in C_{\eta}(\mathbb{R},\mathbb{C}^{m})$, $x\in\mathbb{R}$,
and estimate

(34) $\displaystyle|(L\varphi)(x)|_{\infty}$
$\displaystyle\leq\max_{i\in\mathbb{N}_{m}}\sum_{j\in\mathbb{N}_{m}}\Big{(}|\varphi_{i}(x)|+|\varphi_{j}(x)|\Big{)}$
$\displaystyle+\max_{i\in\mathbb{N}_{m}}\sum_{j\in\mathbb{N}_{m}}\int_{cT_{ji}}^{\infty}e^{-y/c}|w(y)\psi_{ij}(y)|\;|\varphi_{i}(x)-\varphi_{j}(x-y)|\,dy.$

For the first summands in Eq. 34 we find

(35) $|\varphi_{i}(x)|+|\varphi_{j}(x)|\leq 2|\varphi(x)|_{\infty}\leq
2e^{\eta|x|}\|\varphi\|_{C_{m,\eta}}$

For the second summands in Eq. 34, we estimate

(36) $\displaystyle\int_{cT_{ji}}^{\infty}e^{-y/c}$
$\displaystyle|w(y)\psi_{ij}(y)|\;|\varphi_{i}(x)-\varphi_{j}(x-y)|\,dy$ (by
Eq. 33) $\displaystyle\leq
K\int_{cT_{ji}}^{\infty}|w(y)|\,|\varphi_{i}(x)-\varphi_{j}(x-y)|\,dy$
$\displaystyle\leq
K\|\varphi\|_{C_{m,\eta}}\int_{cT_{ji}}^{\infty}|w(y)|\big{(}e^{\eta|x|}+e^{\eta|x-y|}\big{)}\,dy$
$\displaystyle\leq
Ke^{\eta|x|}\|\varphi\|_{C_{m,\eta}}\int_{cT_{ji}}^{\infty}|w(y)|\big{(}1+e^{\eta|y|}\big{)}\,dy$
$\displaystyle\leq
Ke^{\eta|x|}\|\varphi\|_{C_{m,\eta}}\bigg{(}\|w\|_{L^{1}_{\eta}}+\int_{cT_{ji}}^{\infty}|w(y)|e^{\eta|y|}\,dy\bigg{)}$
($w$ even) $\displaystyle\leq
Ke^{\eta|x|}\|\varphi\|_{C_{m,\eta}}\bigg{(}\|w\|_{L^{1}_{\eta}}+2\int_{-\infty}^{\infty}|w(y)|e^{\eta
y}\,dy\bigg{)}$ $\displaystyle\leq
3K\|w\|_{L^{1}_{\eta}}e^{\eta|x|}\|\varphi\|_{C_{m,\eta}}$

Combining Eqs. 34 to 36 we obtain

$\displaystyle\|L\varphi\|_{C_{m,\eta}}$
$\displaystyle=\sup_{x\in\mathbb{R}}e^{-\eta|x|}|L\varphi(x)|_{\infty}$
$\displaystyle\leq m(2+3K\|w\|_{L^{1}_{\eta}})\|\varphi\|_{C_{m,\eta}}$
$\displaystyle:=\kappa_{m,\eta}\|\varphi\|_{C_{m,\eta}},$

which concludes the proof of part 1.

_Part 2_. We set $u=\tau+\varepsilon\varphi\in
C_{\eta}(\mathbb{R},\mathbb{R}^{m})$, for $0<\varepsilon\ll 1$. By main
hypothesis $u$ and $\tau$ are firing functions of two distinct $m$-spike CIFM
solutions. We claim that this implies Equation 26. Indeed, since $u$ is a
firing function, then $V_{m}u=1$ on $\mathbb{F}_{u}$, that is

(37)
$1=I+\sum_{j\in\mathbb{N}_{m}}\Big{(}(Su_{i})(x,u_{j}(x))+(Ru_{j})(x,u_{i}(x))\Big{)},\qquad(i,x)\in\mathbb{N}_{m}\times\mathbb{R}.$

We obtain

$\displaystyle(Su_{j})({\mkern 2.0mu\cdot\mkern 2.0mu},u_{i})$
$\displaystyle=\int_{-\infty}^{\infty}w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\int_{-\infty}^{0}e^{z}\alpha(z+u_{i}-u_{j}(y))\,dz\,dy$
$\displaystyle=(S\tau_{j})({\mkern 2.0mu\cdot\mkern
2.0mu},\tau_{i})+\varepsilon\int_{-\infty}^{\infty}w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\big{(}\varphi_{i}-\varphi_{j}(y)\big{)}\int_{-\infty}^{0}e^{z}\alpha^{\prime}(z+\tau_{i}-\tau_{j}(y))\,dz\,dy$
$\displaystyle+O(\varepsilon^{2}),$

where we have denoted by $\alpha^{\prime}=p^{\prime}H+p\delta$ the
distributional derivative of $\alpha$. We now manipulate the integral in the
previous equation as follows

$\begin{aligned} \int_{-\infty}^{\infty}&w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\big{(}\varphi_{i}-\varphi_{j}(y)\big{)}\int_{-\infty}^{0}e^{z}\alpha^{\prime}(z+\tau_{i}-\tau_{j}(y))\,dz\,dy\\\
&=\int_{-\infty}^{\infty}w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\big{(}\varphi_{i}-\varphi_{j}(y)\big{)}\int_{-\infty}^{\tau_{i}-\tau_{j}(y)}e^{z+\tau_{i}-\tau_{j}(y)}\alpha^{\prime}(z)\,dz\,dy\\\
&=\int_{-\infty}^{\infty}w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\big{(}\varphi_{i}-\varphi_{j}(y)\big{)}\int_{-\infty}^{\tau_{i}-\tau_{j}(y)}e^{z+\tau_{i}-\tau_{j}(y)}\big{(}p^{\prime}(z)H(z)+p(z)\delta(z)\big{)}\,dz\,dy\\\
&=\int_{-\infty}^{\tau^{-1}_{j}(\tau_{i})}w({\mkern 2.0mu\cdot\mkern
2.0mu}-y)\big{(}\varphi_{i}-\varphi_{j}(y)\big{)}e^{\tau_{j}(y)-\tau_{i}}\Big{(}p(0)+\int_{0}^{\tau_{i}-\tau_{j}(y)}e^{z}p^{\prime}(z)\,dz\Big{)}\,dy\\\
\end{aligned},$

hence for all $i,j\in\mathbb{N}_{m}$ we obtain

$(Su_{j})({\mkern 2.0mu\cdot\mkern 2.0mu},u_{i})=(S\tau_{j})({\mkern
2.0mu\cdot\mkern 2.0mu},\tau_{i})+\varepsilon
e^{T_{ji}}\int_{-\infty}^{\infty}e^{-y/c}w(y)\psi_{ij}(y)\big{[}\varphi_{i}-\varphi_{j}({\mkern
2.0mu\cdot\mkern 2.0mu}-y)\big{]}\,dy+O(\varepsilon^{2})$

For the reset operator, we obtain, for all $i,j\in\mathbb{N}_{m}$

$\displaystyle Ru_{j}({\mkern 2.0mu\cdot\mkern 2.0mu},u_{i}^{-})$
$\displaystyle=-\lim_{\kappa\to 0^{+}}\exp(-u_{i}+\kappa+u_{j})H(u_{i}-\kappa-
u_{j})$ $\displaystyle=R\tau_{j}({\mkern 2.0mu\cdot\mkern
2.0mu},\tau_{i}^{-})+\varepsilon\lim_{\kappa\to
0^{+}}\exp(-T_{ij}+\kappa)(\varphi_{i}-\varphi_{j})H(T_{ij})+O(\varepsilon^{2})$
$\displaystyle=R\tau_{j}({\mkern 2.0mu\cdot\mkern
2.0mu},\tau_{i}^{-})+\varepsilon\exp(T_{ji})(\varphi_{i}-\varphi_{j})1_{j<i}+O(\varepsilon^{2}).$

Combining Eq. 37 with the expansions obtained for $S$ and $R$, exploiting the
condition $V_{m}\tau=1$ on $\mathbb{F}_{\tau}$, and dividing by $\varepsilon
e^{T_{ji}}$ we obtain

$0=(L\varphi)_{i}+O(\varepsilon)\qquad\text{on $\mathbb{R}$, for all
$i\in\mathbb{N}_{m}$}.$

which implies Equation 32.

##  Appendix E Proof of Lemma 4.5

######  Lemma E.1 (Restatement of Lemma 4.5).

Assume Section 3.2, let $(c,T)$ be coarse variables of a $\textrm{TW}_{m}$,
and let $\mathbb{D}_{a,b}=\\{z\in\mathbb{C}\colon a\leq\operatorname{Re}z\leq
b\\}$. Further, let $E$ be the complex-valued function

(38) $E\colon\mathbb{D}_{-\eta,\eta}\to\mathbb{C},\quad z\mapsto\det[D-M(z)],$

where $M\in\mathbb{C}^{m\times m}$,
$D=\operatorname{diag}(D_{1},\ldots,D_{m})\in\mathbb{R}^{m\times m}$, are the
matrices with elements

$M_{ij}(z)=e^{T_{ji}}\bigg{[}1_{j<i}+\int_{cT_{ji}}^{\infty}e^{-(z+1/c)y}w(y)\psi_{ij}(y)\,dy\bigg{]},\qquad
D_{i}=\sum_{k\in\mathbb{N}_{m}}M_{ik}(0),$

respectively, then:

  1. 1.

If $\lambda$ is a root of $E$, then its complex conjugate $\lambda^{*}$ is
also a root of $E$, and there exists a nonzero $\Phi\in\ker[D-M(\lambda)]$
such that $\Phi e^{\lambda x},\Phi^{*}e^{\lambda^{*}x}\in\ker{L}$, where $L$
is defined as in Lemma 4.1.

  2. 2.

$E$ has a root at $0$. $\textrm{TW}_{m}$ is linearly stable (but not
asymptotically linearly stable) to perturbations $\varphi\colon x\mapsto\kappa
v$, where $\kappa\in\mathbb{R}\setminus\\{0\\}$ and
$v=(1,\ldots,1)\in\mathbb{R}^{m}$.

  3. 3.

If $\lambda$ is a root of $E$ in
$\mathbb{D}_{-\eta,0}\setminus\mathrm{i}\mkern 1.0mu\mathbb{R}$, then
$\textrm{TW}_{m}$ is linearly asymptotically stable to perturbations $\Phi
e^{\lambda x}+\Phi^{*}e^{\lambda^{*}x}$.

######  Proof E.2.

_Part 1_. We observe that $D$ has purely real entries, and a direct
calculation shows $M(z^{*})=M^{*}(z)$. If $E(\lambda)=0$, then there exists
$\Phi\in\mathbb{C}^{m}\setminus\\{0\\}$ such that $D\Phi=M(\lambda)\Phi$, that
is, $\Phi\in\ker[D-M(\lambda)]$. Taking the complex conjugate we obtain
$D\Phi^{*}=M^{*}(\lambda)\Phi^{*}=M(\lambda^{*})\Phi^{*}$, hence
$E(\lambda^{*})=0$, therefore $\lambda^{*}$ is also a root.

We now set $\varphi=\Phi e^{\lambda x}$ which is in
$C_{\eta}(\mathbb{R},\mathbb{C}^{m})$ because
$\lambda\in\mathbb{D}_{-\eta,\eta}$, and we obtain

$(L\varphi)(x)=e^{\lambda x}[D-M(\lambda)]\Phi=0\qquad x\in\mathbb{R}$

because $\Phi\in\ker[D-M(\lambda)]$. The previous identity implies $\Phi
e^{\lambda x},\Phi^{*}e^{\lambda^{*}x}\in\ker{L}$.

_Part 2_. By definition of $D$ and $M$ we have $[D-M(0)]v=0$, hence $v$ is in
the kernel of $D-M(0)$ and $E(0)=0$. We fix
$\kappa\in\mathbb{R}\setminus\\{0\\}$, use part 1 with $\lambda=0$,
$\Phi=\kappa v$, and deduce that the mapping $\varphi\colon x\mapsto\kappa v$,
which is an element of $C_{\eta}(\mathbb{R},\mathbb{R}^{m})$, is in $\ker L$.
Since $\varphi_{i}(x)\equiv\kappa$, for all $i\in\mathbb{N}_{m}$,
$\textrm{TW}_{m}$ is linearly stable according to Definition 4.4. However,
$|\varphi(x)|_{\infty}\to\kappa\neq 0$ as $x\to\infty$, so $\textrm{TW}_{m}$
is not asymptotically linearly stable.

_Part 3_. Let $\lambda=\mu+\mathrm{i}\mkern 1.0mu\omega$. By main hypothesis
$\mu<0$. From part 1 we deduce that there exists $\Phi\in\ker[D-M(\lambda)]$,
such that $\varphi(x)=\Phi e^{\lambda x}+\Phi^{*}e^{\lambda^{*}x}\in\ker L$.
We note that $|\Phi|_{\infty}$ can be fixed to an arbitrary nonzero constant,
and we bound $\varphi_{i}$ as follows

(39) $|\varphi_{i}(x)|\leq 2|\Phi|_{\infty}e^{\mu x}\leq
K|\Phi|_{\infty},\quad
K=2\max_{j\in\mathbb{N}_{m}}e^{-cT_{j}\mu},\quad(i,x)\in\mathbb{N}_{m}\times[-cT_{i},\infty).$

We now fix $\varepsilon>0$. The bound Eq. 39 and the choices
$\delta=\varepsilon$ and $|\Phi|_{\infty}<\varepsilon/K$ imply that
$\textrm{TW}_{m}$ is linearly stable, according to Definition 4.4. Using again
Eq. 39 we obtain

$\lim_{x\to\infty}|\varphi(x)|_{\infty}\leq
2|\Phi|_{\infty}\lim_{x\to\infty}e^{\mu x}=0$

therefore $\textrm{TW}_{m}$ is linearly asymptotically stable.

##  Appendix F Two-parameter continuation of bifurcations

Grazing points are found generically as a secondary control parameter, say
$\gamma$, is varied. It is possible to perform a $2$-parameter continuation of
the grazing point in the $(\beta,\gamma)$-plane by continuing in $\gamma$
solutions the following problem:

######  Problem F.1 (Grazing point computation).

Find
$(c,T_{1},\ldots,T_{m},T_{G},\beta_{G})\in\mathbb{R}_{>0}\times\mathbb{R}^{m+2}$
such that $T_{1}<\cdots<T_{m}<T_{G}$ and

(40) $\displaystyle T_{1}=0,$ (41)
$\displaystyle\nu([cT_{i}]_{-};c,T_{1},\ldots,T_{m},\beta_{G})=1,\quad
i\in\mathbb{N}_{m},$ (42)
$\displaystyle\nu(cT_{G};c,T_{1},\ldots,T_{m},\beta_{G})=1,$ (43)
$\displaystyle\nu^{\prime}(cT_{G};c,T_{1},\ldots,T_{m},\beta_{G})=0.$

We note that we have exposed the dependence of $\nu$ and $\nu^{\prime}$ on
$\beta$ in the previous problem. Unlike the numerical continuation presented
in the main text, here $\beta$ is a free parameter, which is determined also
by Newton’s method.

Similarly, we can trace loci of Hopf bifurcations in the
$(\gamma,\beta)$-plane, by solving

######  Problem F.2 (Computation of Hopf bifurcations).

Find
$(c,T_{1},\ldots,T_{m},\beta_{\textrm{HB}},\omega_{\textrm{HB}})\in\mathbb{R}_{>0}
\times\mathbb{R}^{m+2}$ such that $T_{1}<\cdots<T_{m}$ and

(44) $\displaystyle T_{1}=0,$ (45)
$\displaystyle\nu([cT_{i}]_{-};c,T_{1},\ldots,T_{m},\beta_{\textrm{HB}})=1,\quad
i\in\mathbb{N}_{m},$ (46) $\displaystyle
E(i\omega_{\textrm{HB}};c,T_{1},\ldots,T_{m},\beta_{\textrm{HB}})=0.$

Figure 11: Bifurcation diagram showing wave speed branches $c$ for
$\textrm{TW}_{m}$ states with $m=1,2,4,8,16$, with purely excitatory kernel,
$a_{1}=2$, $b_{1}=5$, $a_{2}=0$. (a) External current is fixed at $I=0.82$ and
solutions continued in $\beta$, here $m=1,2,4,8$ branches switch stability at
a fold, but for $m=16$ the branch terminates at a grazing point as seen in the
lateral-inhibition case previously. (b) The synaptic processing rate is fixed
at $\beta=4$ and the waves continued in $I$. Note that in both cases a greater
$m$ leads to larger speeds. (c) Synaptic profiles for $\textrm{TW}_{8}$
solutions labelled (1–5) in (a). Smaller $\beta$ values result in more
localised profiles.

## References

  * [1] S.-i. Amari, Dynamics of pattern formation in lateral-inhibition type neural fields, Biolological Cybernetics, 27 (1977), pp. 77–87.
  * [2] D. Avitabile, Numerical computation of coherent structures in spatially-extended systems, May 2020, https://doi.org/10.5281/zenodo.3821169, https://doi.org/10.5281/zenodo.3821169.
  * [3] D. Avitabile and K. C. A. Wedgwood, Macroscopic coherent structures in a stochastic neural network: from interface dynamics to coarse-grained bifurcation analysis, Journal of Mathematical Biology, (2017), pp. 1–44, https://doi.org/10.1007/s00285-016-1070-9.
  * [4] J. Baladron, D. Fasoli, O. Faugeras, and J. Touboul, Mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons, The Journal of Mathematical Neuroscience, 2 (2012), p. 10.
  * [5] D. Barkley, Simplifying the complexity of pipe flow, Physical Review E, 84 (2011), p. 016309.
  * [6] D. Barkley, Pipe flow as an excitable medium, Revista Cubana de Física, 29 (2012), pp. 1–27.
  * [7] D. Barkley, Theoretical perspective on the route to turbulence in a pipe, Journal of Fluid Mechanics, 803 (2016).
  * [8] D. Barkley, B. Song, V. Mukund, G. Lemoult, M. Avila, and B. Hof, The rise of fully turbulent flow, Nature, 526 (2015), pp. 550–553.
  * [9] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, Understanding the dynamics of biological and neural oscillator networks through mean-field reductions: a review, Journal of Mathematical Neuroscience, 10 (2020), pp. 1–43.
  * [10] P. C. Bressloff, Mean-field theory of globally coupled integrate-and-fire neural oscillators with dynamic synapses, Physical Review E, 60 (1999), p. 2160.
  * [11] P. C. Bressloff, Synaptically generated wave propagation in excitable neural media, Physical Review Letters, 82 (1999), pp. 2979–2982.
  * [12] P. C. Bressloff, Traveling waves and pulses in a one-dimensional network of excitable integrate-and-fire neurons, Journal of Mathematical Biology, 40 (2000), pp. 169–198.
  * [13] P. C. Bressloff, Waves in neural media, Lecture Notes on Mathematical Modelling in the Life Sciences, Springer, New York, (2014).
  * [14] P. C. Bressloff and S. Coombes, Dynamics of strongly-coupled spiking neurons., Neural Computation, 12 (2000), pp. 91–129.
  * [15] A. N. Burkitt, A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input, Biological Cybernetics, 95 (2006), pp. 1–19, https://doi.org/10.1007/s00422-006-0068-6.
  * [16] Á. Byrne, D. Avitabile, and S. Coombes, Next-generation neural field model: the evolution of synchrony within patterns and waves, Physical Review E, 99 (2019), p. 012313.
  * [17] M. Camperi and X.-J. Wang, A model of visuospatial working memory in prefrontal cortex: recurrent network and cellular bistability, Journal of Computational Neuroscience, 5 (1998), pp. 383–405.
  * [18] A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang, Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model, Cerebral Cortex, 10 (2000), pp. 910–923.
  * [19] A. Compte, M. V. Sanchez-Vives, D. A. McCormick, and X.-J. Wang, Cellular and network mechanisms of slow oscillatory activity (¡ 1 Hz) and wave propagations in a cortical network model, Journal of Neurophysiology, 89 (2003), pp. 2707–2725.
  * [20] S. Coombes, P. beim Graben, and R. Potthast, Tutorial on Neural Field Theory, in Neural Fields, Springer, Berlin, Heidelberg, Berlin, Heidelberg, 2014, pp. 1–43.
  * [21] S. Coombes, P. beim Graben, R. Potthast, and J. Wright, Neural fields: Theory and Applications, Springer, 2014.
  * [22] S. Coombes and P. C. Bressloff, Saltatory waves in the spike-diffuse-spike model of active dendritic spines, Physical Review Letters, 91 (2003), pp. 81–4.
  * [23] S. Coombes and C. Laing, Pulsating fronts in periodically modulated neural field models, Physical Review E, 83 (2011), p. 011912.
  * [24] S. Coombes, H. Schmidt, and D. Avitabile, Spots: breathing, drifting and scattering in a neural field model, in Neural fields, Springer, Heidelberg, Berlin, Heidelberg, 2014, pp. 187–211.
  * [25] A. Darbyshire and T. Mullin, Transition to turbulence in constant-mass-flux pipe flow, Journal of Fluid Mechanics, 289 (1995), pp. 83–114.
  * [26] F. Delarue, J. Inglis, S. Rubenthaler, E. Tanré, et al., Global solvability of a networked integrate-and-fire model of mckean–vlasov type, The Annals of Applied Probability, 25 (2015), pp. 2096–2133.
  * [27] M. di Bernardo, C. J. Budd, A. R. Champneys, P. Kowalczyk, A. B. Nordmark, G. O. Tost, and P. T. Piiroinen, Bifurcations in Nonsmooth Dynamical Systems, SIAM Review, 50 (2008), pp. 629–701.
  * [28] J. R. Dormand and P. J. Prince, A family of embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics, 6 (1980), pp. 19–26, https://doi.org/10.1016/0771-050x(80)90013-3.
  * [29] B. Eckhardt, H. Faisst, S. A., and J. Schumacher, Turbulence transition in shear flows., in Advances in Turbulence IX: Proc. Ninth European Turbulence Conference, Barcelona, I. P. Castro, P. E. Hancock, and T. G. Thomas, eds., 2002, pp. 701–708.
  * [30] B. Ermentrout, Neural networks as spatio-temporal pattern-forming systems, Reports on Progress in Physics, 61 (1998), p. 353.
  * [31] B. Ermentrout, The analysis of synaptically generated traveling waves, Journal of Computational Neuroscience, 5 (1998), pp. 191–208.
  * [32] B. Ermentrout, J. Rubin, and R. Oşan, Regular traveling waves in a one-dimensional network of theta neurons, SIAM Journal on Applied Mathematics, 62 (2002), pp. 1197–1221.
  * [33] G. B. Ermentrout, S. E. Folias, and Z. P. Kilpatrick, Spatiotemporal pattern formation in neural fields with linear adaptation, in Neural Fields, Springer, Berlin, Heidelberg, Berlin, Heidelberg, 2014, pp. 119–151.
  * [34] G. B. Ermentrout and N. Kopell, Parabolic bursting in an excitable system coupled with a slow oscillation, SIAM Journal on Applied Mathematics, 46 (1986), pp. 233–253.
  * [35] G. B. Ermentrout and D. H. Terman, Mathematical Foundations of Neuroscience, vol. 35, Springer Science & Business Media, 2010.
  * [36] J. M. Esnaola-Acebes, A. Roxin, D. Avitabile, and E. Montbrió, Synchrony-induced modes of oscillation of a neural field model, Physical Review E, 96 (2017), p. 052407.
  * [37] H. Faisst and B. Eckhardt, Traveling waves in pipe flow, Physical Review Letters, 91 (2003), p. 224502.
  * [38] S. E. Folias and P. C. Bressloff, Breathing pulses in an excitatory neural network, SIAM Journal on Applied Dynamical Systems, 3 (2004), pp. 378–407.
  * [39] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, Neuronal dynamics: from single neurons to networks and models of cognition, Cambridge University Press, 2014.
  * [40] W. Gerstner, J. L. van Hemmen, and J. D. Cowan, What matters in neuronal locking?, dx.doi.org, 8 (2008), pp. 1653–1676.
  * [41] J. F. Gibson, J. Halcrow, and P. Cvitanović, Equilibrium and travelling-wave solutions of plane couette flow, Journal of Fluid Mechanics, 638 (2009), pp. 243–266.
  * [42] D. Golomb and G. B. Ermentrout, Continuous and lurching traveling pulses in neuronal networks with delay and spatially decaying connectivity, Proceedings of the National Academy of Sciences, 96 (1999), pp. 13480–13485.
  * [43] L. R. González-Ramírez, O. J. Ahmed, S. S. Cash, C. E. Wayne, and M. A. Kramer, A biologically constrained, mathematical model of cortical wave propagation preceding seizure termination, PLoS Computational Biology, 11 (2015), pp. e1004065–34.
  * [44] A. Granados, L. Alseda, and M. Krupa, The period adding and incrementing bifurcations: from rotation theory to applications, SIAM Review, 59 (2017), pp. 225–292.
  * [45] B. Hof, A. Juel, and T. Mullin, Scaling of the turbulence transition threshold in a pipe, Physical Review Letters, 91 (2003), p. 244502.
  * [46] E. Hopf, A mathematical example displaying features of turbulence, Communications on Pure and Applied Mathematics, 1 (1948), pp. 303–322.
  * [47] X. Huang, W. C. Troy, Q. Yang, H. Ma, C. R. Laing, S. J. Schiff, and J.-Y. Wu, Spiral waves in disinhibited mammalian neocortex, J. Neurosci., 24 (2004), pp. 9897–9902.
  * [48] J. Inglis and J. MacLaurin, A general framework for stochastic traveling waves and patterns, with application to neural field equations, SIAM Journal on Applied Dynamical Systems, 15 (2016), pp. 195–234.
  * [49] E. M. Izhikevich, Dynamical systems in neuroscience, MIT press, 2007\.
  * [50] Z. P. Kilpatrick and B. Ermentrout, Wandering bumps in stochastic neural fields, SIAM Journal on Applied Dynamical Systems, 12 (2013), pp. 61–94.
  * [51] S. S. Kim, H. Rouault, S. Druckmann, and V. Jayaraman, Ring attractor dynamics in the Drosophila central brain, Science, 356 (2017), pp. 849–853.
  * [52] J. J. Knierim and K. Zhang, Attractor dynamics of spatially correlated neural activity in the limbic system, Annual Review of Neuroscience, 35 (2012), pp. 267–285.
  * [53] C. R. Laing, Exact neural fields incorporating gap junctions, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 1899–1929.
  * [54] C. R. Laing and C. C. Chow, Stationary Bumps in Networks of Spiking Neurons, Neural Computation, 13 (2001), pp. 1473–1494.
  * [55] C. R. Laing and O. Omel’chenko, Moving bumps in theta neuron networks, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30 (2020), p. 043117.
  * [56] L. D. Landau, On the problem of turbulence, in Doklady Akademii Nauk USSR, vol. 44, 1944, p. 311.
  * [57] L. Lapicque and M. Lapicque, Recherches quantitatives sur l’excitation électrique des nerfs traitée comme une polarisastion, Journal de Physiologie et de Pathologie Général, 9 (1907), pp. 620–635.
  * [58] T. B. Luke, E. Barreto, and P. So, Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons, Neural Computation, 25 (2013), pp. 3207–3234.
  * [59] J. MacLaurin and P. Robinson, Determination of effective brain connectivity from activity correlations, Physical Review E, 99 (2019), p. 042404.
  * [60] P. Manneville, On the transition to turbulence of wall-bounded flows in general, and plane couette flow in particular, European Journal of Mechanics-B/Fluids, 49 (2015), pp. 345–362.
  * [61] A. Meseguer and L. N. Trefethen, Linearized pipe flow to Reynolds number $10^{7}$, Journal of Computational Physics, 186 (2003), pp. 178–197.
  * [62] R. E. Mirollo and S. H. Strogatz, Synchronization of pulse-coupled biological oscillators, SIAM Journal on Applied Mathematics, 50 (2006), pp. 1645–1662.
  * [63] E. Montbrió, D. Pazó, and A. Roxin, Macroscopic description for networks of spiking neurons, Physical Review X, 5 (2015), p. 021028.
  * [64] R. Oşan, R. Curtu, J. Rubin, and B. Ermentrout, Multiple-spike waves in a one-dimensional integrate-and-fire neural network, Journal of Mathematical Biology, 48 (2004), pp. 243–274.
  * [65] R. Oşan and B. Ermentrout, The evolution of synaptically generated waves in one- and two-dimensional domains, Physica D: Nonlinear Phenomena, 163 (2002), pp. 217–235.
  * [66] C. C. Pringle, Y. Duguet, and R. R. Kerswell, Highly symmetric travelling waves in pipe flow, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367 (2009), pp. 457–472.
  * [67] J. Rankin, D. Avitabile, J. Baladron, G. Faye, and D. J. Lloyd, Continuation of localized coherent structures in nonlocal neural field equations, SIAM Journal on Scientific Computing, 36 (2014), pp. B70–B93.
  * [68] A. D. Redish, A. N. Elga, and D. S. Touretzky, A coupled attractor model of the rodent head direction system, Network: Computation in Neural Systems, 7 (1996), pp. 671–685.
  * [69] O. Reynolds, An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels, Philosophical Transactions of the Royal Society of London, (1883), pp. 935–982.
  * [70] K. A. Richardson, S. J. Schiff, and B. J. Gluckman, Control of Traveling Waves in the Mammalian Cortex, Physical Review Letters, 94 (2005), p. 028103.
  * [71] D. Ruelle and F. Takens, On the nature of turbulence, Les rencontres physiciens-mathématiciens de Strasbourg-RCP25, 12 (1971), pp. 1–44.
  * [72] L. Sacerdote and M. T. Giraudo, Stochastic integrate and fire models: a review on mathematical methods and their applications, in Stochastic Biomathematical Models, Springer, 2013, pp. 99–148.
  * [73] H. Salwen, F. W. Cotton, and C. E. Grosch, Linear stability of poiseuille flow in a circular pipe, Journal of Fluid Mechanics, 98 (1980), pp. 273–284.
  * [74] H. Schmidt and D. Avitabile, Bumps and oscillons in networks of spiking neurons, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30 (2020), p. 033133.
  * [75] A. Schmiegel and B. Eckhardt, Fractal stability border in plane couette flow, Physical Review Letters, 79 (1997), p. 5250.
  * [76] L. F. Shampine and M. W. Reichelt, The MATLAB ODE suite, SIAM Journal on Scientific Computing, 18 (1997), pp. 1–22.
  * [77] M. L. Steyn-Ross, D. A. Steyn-Ross, J. W. Sleigh, and D. R. Whiting, Theoretical predictions for spatial covariance of the electroencephalographic signal during the anesthetic-induced phase transition: increased correlation length and emergence of spatial self-organization, Physical Review E, 68 (2003), p. 021902.
  * [78] H. C. Tuckwell, Introduction to Theoretical Neurobiology. Volume 1: Linear Cable Theory and Dendritic Structure, vol. 8, Cambridge University Press, 1988.
  * [79] D. Turner-Evans, S. Wegener, H. Rouault, R. Franconville, T. Wolff, J. D. Seelig, S. Druckmann, and V. Jayaraman, Angular velocity integration in a fly heading circuit, eLife, 6 (2017), p. e23496.
  * [80] C. W. van Doorne and J. Westerweel, The flow structure of a puff, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367 (2009), pp. 489–507.
  * [81] C. van Vreeswijk, Partial synchronization in populations of pulse-coupled oscillators, Physical Review E, 54 (1996), pp. 5522–5537.
  * [82] C. van Vreeswijk, L. F. Abbott, and G. B. Ermentrout, When inhibition not excitation synchronizes neural firing, Journal of Computational Neuroscience, 1 (1994), pp. 313–321.
  * [83] H. Wedin and R. R. Kerswell, Exact coherent structures in pipe flow: Travelling wave solutions, Journal of Fluid Mechanics, 508 (2004), pp. 333–371.
  * [84] H. R. Wilson and J. D. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Kybernetik, 13 (1973), pp. 55–80.
  * [85] K. Wimmer, D. Q. Nykamp, C. Constantinidis, and A. Compte, Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory, Nature Neuroscience, 17 (2014), pp. 431–439.
  * [86] M. Wolfrum, O. E. Omel’chenko, and J. Sieber, Regular and irregular patterns of self-localized excitation in arrays of coupled phase oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science, 25 (2015), pp. 053113–8.
  * [87] K. Zhang, Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory, Journal of Neuroscience, 16 (1996), pp. 2112–2126.