∎

11institutetext: Peter Ashwin 22institutetext: Center for Systems, Dynamics
and Control, Department of Mathematics, University of Exeter, Exeter EX4 4QF,
UK 22email<EMAIL_ADDRESS>33institutetext: Claire Postlethwaite
44institutetext: Department of Mathematics, University of Auckland, Auckland,
1142, New Zealand.

# Excitable Networks for Finite State Computation with Continuous Time
Recurrent Neural Networks

Peter Ashwin  Claire Postlethwaite

(Received: date / Accepted: date)

###### Abstract

Continuous time recurrent neural networks (CTRNN) are systems of coupled
ordinary differential equations that are simple enough to be insightful for
describing learning and computation, from both biological and machine learning
viewpoints. We describe a direct constructive method of realising finite state
input-dependent computations on an arbitrary directed graph. The constructed
system has an excitable network attractor whose dynamics we illustrate with a
number of examples. The resulting CTRNN has intermittent dynamics:
trajectories spend long periods of time close to steady-state, with rapid
transitions between states. Depending on parameters, transitions between
states can either be _excitable_ (inputs or noise needs to exceed a threshold
to induce the transition), or _spontaneous_ (transitions occur without input
or noise). In the excitable case, we show the threshold for excitability can
be made arbitrarily sensitive.

###### Keywords:

Continuous Time Recurrent Neural Network Nonlinear Dynamics Excitable Network
Attractor

††journal: Biological Cybernetics

##  1 Introduction

It is natural to try to understand computational properties of neural systems
through the paradigm of network dynamical systems, where a number of
dynamically simple units (i.e. with attracting equilibria or periodic orbits)
interact to give computation as an emergent property of the system. This is
the as much the case for biological models of information processing, pattern
generation and decision making as it is for artificial neural networks
inspired by these. A variety of specific models have been developed to
describe the dynamics and training of recurrent networks comprised of coupled
neurons in biological and artificial settings. The particular challenge that
we address here is the construction of arbitrarily complex, but specified,
dynamical structures that enable discrete (finite-state) computation in an
input-driven system that is nonetheless continuous in both state and time.
Clearly, invariant objects of the autonomous system (such as equilibria,
periodic orbits and chaotic attractors) only form part of the picture and an
input-dependent (non-autonomous) approach such as manjunath2012theory  is
needed.

To help understand the response of systems to inputs the authors introduced in
Ashwin2016a  a notion of a “network attractor”, namely an invariant object in
phase space that may contain several local invariant sets, but also systems of
interconnections between them. This generalises the notion of “heteroclinic
network” and “winnerless competition/stable heteroclinic channels”
afraimovich2004origin  which have been used to describe a range of sequence
generation, computational and cognitive effects in neural systems: see for
example rabinovich2001dynamical ; rabinovich2006generation ;
afraimovich2004heteroclinic ; rabinovich2020sequential ; hutt2017sequences .
These models connect computational states (represented as saddles in the
dynamics). In the presence of inputs or noise, the switching between saddles
is a useful model for spontaneous computational processes but there are two
problems. One of these is that the states are dynamically unstable and so
there are spontaneous transitions that are not noise driven. The other is that
heteroclinic chains are destroyed by arbitrarily small perturbations, unless
there are special structures in phase space (symmetries or invariant
subspaces).

Network attractors joining stable equilibria Ashwin2016a  overcome both of
these issues. In this paper we show that arbitrarily complex network
attractors exist in idealised class of model is the continuous time recurrent
neural network (CTRNN) beer1995dynamics . A CTRNN is a set of differential
equations each with one scalar variable that represents the level of
activation of a neuron, and feedback via a saturating nonlinearity or
“activation function”. These models have been extensively investigated in the
past decades as simple neurally-inspired systems that can nonetheless (without
input) have complex dynamics by virtue of the nonlinearities present
beer1995dynamics . They can also be trained to perform arbitrarily complex
tasks depending on time-varying input tuci2002evolutionary ;
yamauchi1994sequential . They are frequently used in investigations of
evolutionary robotics blynel2003exploring  and (in various equivalent
formulations chow2019before ) as models for neural behaviour in biological or
cognitive systems. For example, the classical work of Hopfield and Tank
hopfield1985neural  considers such systems with symmetric weights to solve
optimization problems, while more recently, bhowmik2016reservoir  discusses
CTRNN models for episodic memory, and several other biological and cognitive
applications are discussed in nikiforou2019dynamics .

CTRNNs are often referred to as “universal dynamical approximators”
funahashi1993approximation , meaning that the trajectory of a CTRNN can
approximate, to arbitrary precision, any other prescribed (smooth) trajectory
in $\mathbb{R}^{n}$. However, this does not mean that the dynamics of CTRNNs
are simple to understand, or that it is easy to form the above approximation.
It also raises the question of how a “trained” CTRNN performs a complex task.
Gradient descent or more general evolutionary training algorithms train the
network by navigating through a high dimensional landscape of possible
feedback weightings and moving these towards a setting that is sufficiently
optimal for the task. It may be possible to give a clear description of the
resulting nonlinear dynamics of the autonomous (constant input) CTRNN, but we
want to understand not only this but also how inputs affect the state of the
system.

The main theoretical result in this paper focusses on “excitable network
attractors”: these consist of a finite set of local attracting equilibria and
excitable connections between them (see Appendix A). It was demonstrated in
Ashwin2018sensitive  that an excitable network attractor can be used to embed
an arbitrary Turing machine within a class of purpose-designed coupled
dynamical system with two different cell types. Rather than relying on an
optimization approach to design the system, that paper gave a _constructive_
method for designing a realisation of any desired network attractor. However,
this construction required specialist dynamical cells with quite complex
nonlinear couplings between them, and a comparatively large number of cells.
It was recently shown in ceni2020interpreting  that trained RNNs in the form
of echo-state networks can realise excitable networks for certain tasks and
that structural errors in the trained network can explain errors in
imperfectly trained systems.

The current paper demonstrates that CTRNN dynamics is sufficiently rich to
realise excitable networks with arbitrary graph topology, simply by specifying
appropriate connection weights. The construction algorithm in the proof of
Theorem 2.1, assigns one of only four values to each of the connection weights
to realise an arbitrary graph on $N$ vertices as an excitable network on $N$
states (subject to some minor constraints on its connectivity), using a CTRNN
with $N$ cells. The CTRNN we consider in this paper (see for example
beer1995dynamics , which corresponds to a continuous time Hopfield model
hopfield1985neural  in the symmetric coupling case $w_{ij}=w_{ji}$) are
ordinary differential equations

$\dot{y}_{i}=-y_{i}+\sum_{j}w_{ij}\phi(y_{j})+I_{i}(t),$ (1)

where $\mathbf{y}=(y_{1},\ldots y_{N})\in\mathbb{R}^{N}$ is the internal state
of the $N$ cells of the system, $w_{ij}$ is a matrix of connection weights,
$\phi$ is a (sigmoid) activation function that introduces a saturating
nonlinearity into the system and $I_{i}(t)$ is an input. We say system (1) is
input-free if $I_{i}(t)=0$ for all $i$ and $t$.

We consider two cases for $\phi$, a smooth function

$\phi(y)=\phi_{S}(y):=\left[1+\exp\left(-\frac{(y-\theta)}{\epsilon}\right)\right]^{-1},$
(2)

and a piecewise affine function

$\phi(y)=\phi_{P}(y):=\begin{cases}0&y-\theta<-2\epsilon,\\\
(y-\theta)/(4\epsilon)+1/2&|y-\theta|\leq 2\epsilon,\\\
1&y-\theta>2\epsilon,\end{cases}$ (3)

In both cases, $\phi$ is monotonic increasing with

$\lim_{y\rightarrow\infty}\phi(y)=1,~{}\lim_{y\rightarrow-\infty}\phi(y)=0,~{}\phi(\theta)=1/2,$

and a maximum derivative at $y=\theta$, equal to $\frac{1}{4\epsilon}.$ In
both cases, $\epsilon$ and $\theta$ are parameters, and we are interested in
the case $0<\epsilon\ll 1$. In general, the function $\phi$ need not be the
same in every component of (1), but here we make a simplifying assumption that
it is.

Note that both activation functions (2) and (3) have piecewise constant limits
in the singular limit $\epsilon\rightarrow 0$. Such limiting systems are of
Fillipov type and have been explored in various biological contexts,
especially for gene regulatory dynamics. These can also have rich dynamics as
discussed in the literature (for example gouze2003aclass ;
harris2015bifurcations ) but we do not consider this limit here.

The main contribution of Section 2 is to give, in Theorem 2.1, a construction
of a connection weight matrix $w_{ij}$ such that the dynamics of the input-
free system (1) contains (or _realises_ : definition given below) an excitable
network attractor, as defined in Ashwin2016a . We prove this (with details in
Appendix B) for the case of the piecewise affine function $\phi_{P}$. In
section 3 we present evidence that this is also true for the smooth case
$\phi_{S}$ for an open set of parameters. Qualitatively, this means the system
will contain a number of stable equilibrium states, and small inputs (either
deterministic, or noisy) will push the trajectory from one stable equilibrium
into the basin of attraction of another. In this way, transitions can be made
around the network, and the transition time between states tends to be much
smaller than the residence times of the trajectory in neighbourhoods of the
states. In particular, we can choose $w_{ij}$ so that the network attractor
has (almost) any desired topology. In Appendix A we recall formal definitions
of network attractors from Ashwin2016a ; Ashwin2018sensitive .

In section 3 we consider several examples of simple graphs, and demonstrate
that the desired networks do indeed exist in the systems as designed. We also
perform numerical bifurcation analysis to demonstrate the connection between
periodic orbits in the input-free deterministic system (1) and excitable
networks in the same system with additive noise, that is, the system of
stochastic differential equations (SDEs):

$dy_{i}=\left(-y_{i}+\sum_{j}w_{ij}\phi(y_{j})\right)dt+\sigma dW_{i}(t),$ (4)

where $W_{i}(t)$ are independent standard Wiener processes. Here, the noise
plays the role of inputs that propel the trajectory around the network,
although of course this occurs in a random manner. In sections 3.3 and 3.4 we
consider graphs that have multiple edges leading out from a single vertex, and
show that additional equilibria may appear in the network attractor where two
or more cells are active simultaneously. We further show that the existence of
these additional equilibria can be suppressed by choosing one of the
parameters used in the construction of the weight matrix $w_{ij}$ to be
sufficiently large.

Section 4 concludes by relating our results to other notions of sequential
computation. We also conjecture some extensions of the results shown in this
paper.

##  2 Construction of a CTRNN with a network attractor

Let $G$ be an arbitrary directed graph between $N$ vertices, and let $a_{ij}$
be the adjacency matrix of $G$. That is, $a_{ij}=1$ if there is a directed
edge from vertex $i$ to vertex $j$, and $a_{ij}=0$ otherwise.

Let $\Sigma$ be an invariant set for a system of ordinary differential
equations. We say $\Sigma$ is an excitable network that _realises_ a graph $G$
for some amplitude $\delta>0$ if for each vertex $v_{i}$ in $G$ there is a
unique stable equilibrium $\xi_{i}$ in $\Sigma$, and if there is an excitable
connection with amplitude $\delta$ in $\Sigma$ from $\xi_{i}$ to $\xi_{j}$
whenever there is an edge in $G$ from $v_{i}$ to $v_{j}$. The existence of an
excitable connection means that there exists a trajectory with initial
condition within a distance $\delta$ of $\xi_{i}$, which asymptotes in forward
time to $\xi_{j}$ (formal definitions are given in Appendix A).

For the purposes of our construction of a network attractor, we assume that
$G$ contains no loops of order one, no loops of order two, and no
$\Delta$-cliques. Figure 1 shows each of these graph components schematically.

Figure 1:  From left to right, the figures show an order-one loop, an order-
two loop, and a $\Delta$-clique in a directed graph. The construction we
provide realises an arbitrary graph $G$, as long of none of these subgraphs
are present in $G$.

In terms of the adjacency matrix, for $G$ to contain no loops of order one
requires that

$a_{ii}=0\mbox{ for all }i;$ (5)

for $G$ to contain no loops of order two requires that

$a_{ij}a_{ji}=0\mbox{ for all }i,j;$ (6)

and for $G$ to contain no $\Delta$-cliques requires that

$a_{ij}a_{ik}a_{jk}=0\mbox{ for all }i,j,k.$ (7)

In our earlier work, we have demonstrated a network design which can admit
order-two loops and $\Delta$-cliques Ashwin2016a , although this previous
construction is not motivated by neural networks _per se_ , and requires a
higher dimensional system of ODEs (for a given graph) than the one presented
here.

Before we move into the details of the construction, we briefly discuss our
terminology. A graph $G$ has _vertices_ , which correspond to (stable)
_equilibria_ in the phase space of the dynamical system (1). Also within the
phase space, there exist _excitable connections_ (sometimes abbreviated to
_connections_) between the equilibria, which correspond to the _edges_ of the
graph. When a trajectory in the phase space moves between neighbourhoods of
equilibria along a path close to one of these connections, we say that a
_transition_ between the equilibria has occurred. We refer to each of the
components of the dynamical system (1) as a _cell_ , and say that a cell $j$
is _active_ if $\phi(y_{j})$ is close to one.

###  2.1 Realization of arbitrary directed graphs as network attractors

We construct a weight matrix $w_{ij}$ that depends only on the adjacency
matrix $a_{ij}$, and on four parameters $w_{t}$, $w_{s}$, $w_{m}$ and $w_{p}$.
It is given by

$w_{ij}=w_{t}+(w_{s}-w_{t})\delta_{ij}+(w_{p}-w_{t})a_{ji}+(w_{m}-w_{t})a_{ij},$
(8)

where $\delta_{ij}$ is the Kronecker $\delta$. This choice of $w_{ij}$ ensures
that $w_{ii}=w_{s}$, $w_{ij}=w_{p}$ if there is a directed edge in $G$ from
vertex $j$ to $i$ (i.e. $a_{ji}=1$), $w_{ij}=w_{m}$ if there is a directed
edge in $G$ from vertex $i$ to vertex $j$ (i.e. $a_{ij}=1$), and
$w_{ij}=w_{t}$ otherwise. In later sections we allow for different weights
along different edges by allowing $w_{p}$ to depend on $i$ and $j$ (i.e.
$w_{p}=w_{p}^{ij}$). We give an overview of how each of the parameters affect
the dynamics of the system in section 2.2.

We write
$\mathbf{w}=(\epsilon,\theta,w_{s},w_{m},w_{t},w_{p})\in\mathbb{R}^{6}$ to be
a vector of all parameter values. The next result shows that for the piecewise
affine activation function and suitable choice of parameters, there is an
embedding of $G$ as an excitable network attractor for the input-free system.

###### Theorem 2.1

For any directed graph $G$ with $N$ vertices containing no loops of order one,
no loops of order two, and no $\Delta$-cliques, and any small enough
$\delta>0$, there is an open set $W_{\mathrm{ex}}\subset\mathbb{R}^{6}$ such
that if the parameters $\mathbf{w}\in W_{\mathrm{ex}}$ then the dynamics of
input-free equation (1) on $N$ cells with piecewise affine activation (3) and
$w_{ij}$ defined by (8) contains an excitable network attractor with threshold
$\delta$, that realises the graph $G$.

Recall that by _realises_ we mean that all edges in the graph are present as
transitions between stable equilibria using perturbations of size at most
$\delta$.

###### Proof

We give the main ideas behind the proof here, deferring some of the details to
Appendix B. We construct an excitable network attractor in $\mathbb{R}^{N}$
for (1) with piecewise activation function (3) and weight matrix (8). For any
$\frac{1}{2}>\delta>0$, we show there exist parameters $\bf{w}$ (with
$\epsilon>0$ small) and stable equilibria $\xi_{k}$ ($k=1,\ldots,N$) that are
connected according to the adjacency matrix $a_{ij}$ by excitable connections
with amplitude $\delta$. Below, we provide an explicit set of parameters that
make such a realisation, and note that the realisation will hold for an open
set of nearby parameters.

We show in Appendix B that the equilibria $\xi_{k}$ have components (cells)
that are close to one of four values $Y_{T}$, $Y_{D}$, $Y_{L}$, $Y_{A}$
related to the edges attached to the corresponding vertex $k$ in the graph
$G$. For any $0<\delta<\frac{1}{2}$, we use the following parameters:

$\epsilon=\dfrac{\delta}{8},~{}\theta=\frac{1}{2},~{}w_{s}=1,~{}w_{t}=0,$ (9)

and then $w_{p}$ and $w_{m}$ are given by

$w_{p}=\theta-\dfrac{\delta}{2},\quad
w_{m}=-(w_{s}-\theta)-\dfrac{\delta}{2}.$ (10)

We then set

$\begin{split}Y_{A}:=w_{s},\ Y_{L}:=w_{p}=\theta-\frac{\delta}{2},\ \\\
Y_{T}:=w_{m}=-(w_{s}-\theta)-\frac{\delta}{2},\ Y_{D}:=w_{t}\end{split}$ (11)

We use square brackets and subscripts to identify the components of points in
phase space, that is, $[\xi_{k}]_{j}$ is the $j$th component of $\xi_{k}$.
Each $\xi_{k}$ has:

  * •

Exactly one cell that is Active: $[\xi_{k}]_{k}=Y_{A}$

  * •

A number of cells that are Leading: $[\xi_{k}]_{j}=Y_{L}$ (if $a_{kj}=1$)

  * •

A number of cells that are Trailing: $[\xi_{k}]_{j}=Y_{T}$ (if $a_{jk}=1$)

  * •

All remaining cells are Disconnected: $[\xi_{k}]_{j}=Y_{D}$
($a_{kj}=a_{jk}=0$).

Note that the requirement of no loops of order one or two and no
$\Delta$-cliques implies that this labelling is well defined.

From equilibrium $\xi_{k}$, there is an excitable connection to any of the
equilibria $\xi_{l}$ with $a_{kl}=1$, that is, any of the Leading cells can
become the Active cell. During a transition, the remaining cells can be
classified into six types, which are identified in figure 2, and depend (for
each $j$) on the values of the four entries in the adjacency matrix $a_{jk}$,
$a_{kj}$, $a_{jl}$ and $a_{lj}$. We label cell $k$ as AT (Active–Trailing) and
cell $l$ as LA (Leading–Active). Note that the lack of two cycles means that
the cases with $a_{jk}=a_{kj}=1$ or $a_{jl}=a_{lj}=1$ (a total of seven
possibilities) cannot occur, and the lack of $\Delta$-cliques mean that the
cases with $a_{jk}=a_{jl}=1$, $a_{kj}=a_{jk}=1$, or $a_{kj}=a_{lj}=1$ also
cannot occur (which includes the cases where a cell would switch from Leading
to Trailing). The remaining six possibilities are listed below.

Figure 2: Schematic diagram showing a transition in the network attractor of
Theorem 2.1. Left: a schematic of the directed graph showing the eight
distinct types of vertex. Right: schematic time series; as the Active cell
changes there will be a transition in the connected Leading and Trailing cells
as shown. Which cell becomes active is controlled by a small perturbation to
the set of currently Leading cells.

  * •

Type DD: $a_{jk}=a_{kj}=a_{jl}=a_{lj}=0$; the cell is Disconnected throughout.

  * •

Type TD: $a_{jk}=1$, $a_{kj}=a_{jl}=a_{lj}=0$; the cell switches from Trailing
to Disconnected.

  * •

Type LD: $a_{kj}=1$, $a_{jk}=a_{jl}=a_{lj}=0$; the cell switches from Leading
to Disconnected.

  * •

Type TL: $a_{jk}=a_{lj}=1$, $a_{kj}=a_{jl}=0$; the cell switches from Trailing
to Leading.

  * •

Type DT: $a_{jl}=1$, $a_{jk}=a_{kj}=a_{lj}=0$; the cell switches from
Disconnected to Trailing.

  * •

Type DL: $a_{lj}=1$, $a_{jk}=a_{kj}=a_{jl}=0$; the cell switches from
Disconnected to Leading.

The right panel in figure 2 shows how a transition from cell AT active to cell
LA active will occur in a general network.

To prove the existence of an excitable connection giving a realisation we
consider a perturbation from $\xi_{k}$ to the point

$\zeta_{k,l}=\xi_{k}+\delta e_{l},$

where $e_{l}$ is the unit basis vector, and we show in Appendix B that, for
small enough $\delta$, $\zeta_{k,l}$ is in the basin of attraction of
$\xi_{l}$ if $a_{kl}=1$. This means there is an excitable connection from
$\xi_{k}$ to $\xi_{l}$ in this case. ∎

We believe that for small enough $\delta$ and suitable choice of weights, the
realisation of $G$ in Theorem 2.1 can be made almost complete in the sense
analogous to the similar notion for heteroclinic networks ashwin2020almost ,
namely that the set

$\bigcup_{\xi\in E}B_{\delta}(\xi)\setminus\Sigma_{E}$

has zero measure. If a network realization is almost complete then almost all
trajectories starting close to some $\xi_{k}$ will remain close, or will
follow a connection corresponding to the realization. We do not have a proof
of this, though Appendix B.2 shows that for small enough $\delta$,
$\zeta_{k,l}$ is in the basin of attraction of $\xi_{l}$ if and only if
$a_{kl}=1$. This does not preclude the possibility that there exist
perturbations in $B_{\delta}(\xi_{i})$ other than $\zeta_{k,l}$ resulting in
trajectories asymptotic to $\xi_{l}$, or indeed to other attractors. Our
numerical investigations suggest that by choosing $w_{t}$ non-zero and large
enough, any connections to other equilibria can be suppressed. This is
discussed more in sections 3.4 and 4.

###  2.2 Excitable networks for smooth nonlinearities

For small enough $\epsilon$, the smooth activation function (2) can be made
arbitrarily close to the piecewise activation function (3), and so we expect
Theorem 2.1 to also apply in the smooth case, but do not give a proof here.
However, throughout section 3 we use the smooth activation function (2) in our
examples. We use

$\begin{split}&\epsilon=0.05,~{}\theta=0.5,~{}w_{s}=1,\\\
&w_{m}=-0.7,w_{p}=0.3,~{}w_{t}=0\end{split}$ (12)

as our default parameter set (compare this choice of parameters with that
given in equations (9) and (10), with $\delta=0.4$), though there will be an
open set of nearby parameters with analogous behaviour. In section 3 we
provide several examples of using this choice of weight matrix to realise a
graph $G$.

From equation (11) we can see that the parameter choices directly affect the
location of the equilibria $\xi_{k}$ in phase space. As we will see in the
following sections, the parameters also have further effects on the dynamics.
In particular, the relative sizes of the parameters $w_{p}$ and $\theta$
determine whether the dynamics are excitable or spontaneous: essentially, for
$\epsilon$ small enough, $w_{p}$ needs to be smaller than $\theta$ to observe
excitable dynamics. If $w_{p}$ is too large, then the equilibria $\xi_{k}$
cease to exist: periodic orbits can exist instead. The parameter $w_{m}$
controls how fast a Trailing cell decays, and the parameter $w_{t}$ controls
the suppression effects when there is more than one Leading cell. We discuss
the effects of $w_{t}$ in more detail in section 3.3.

For a smooth activation function such as (2) that is invertible on its range,
there is a useful change of coordinates to $J_{i}=\phi(y_{i})$ (a similar
transformation is made in beer1995dynamics ). The input-free equations (1)
then become

$\dot{J}_{i}=\frac{1}{\epsilon}J_{i}(1-J_{i})\left(-\phi^{-1}(J_{i})+\sum_{j}w_{ij}J_{j}\right),$
(13)

where each $J_{i}\in(0,1)$ (which is the domain of the function $\phi^{-1}$),
and

$\phi^{-1}(x)=\theta-\epsilon\ln\left(\frac{1-x}{x}\right).$

Each vertex in the graph $G$ is realised in the phase space of the $J_{k}$
variables by a stable equilibrium with one of the $J_{k}$ close to $1$ and the
remainder close to $0$. With a slight abuse of notation, we still refer to
these equilibria as $\xi_{k}$. As we will see in the examples which follow,
the parameters are chosen such that the dynamics are close to a saddle-node
bifurcation; more precisely such that the system is near a codimension $N$
bifurcation where there are $N$ saddle-nodes of the equilibria $\xi_{k}$. This
bifurcation has a global aspect analogous to a saddle-node on an invariant
circle (SNIC) or saddle-node homoclinic bifurcation in that coupling weights
are such that there are global connections between the saddle nodes that
reflect the network structure.

If parameters are on one side of the $N$ saddle-node bifurcations then for
each $k$ there will exist a pair of equilibria, one of which will be stable
(this is the equilibrium $\xi_{k}$), and a nearby saddle. A small perturbation
can then move the trajectory out of the basin of attraction of the stable
equilibrium and effect a transition to another equilibrium: this is an
_excitable_ connection from $\xi_{k}$.

If parameters are on the other side of the $N$ saddle-node bifurcations then
the flow through the corresponding $N$ region of phase space will be slow but
we will still observe intermittent dynamics as the trajectory passes through
this region (strogatz1994nonlinear, , p99). In this case, we refer to a region
where (a) there is a unique local minimum of $|\dot{\mathbf{y}}|$ and (b) a
large subset of initial conditions in this region pass close to this minimum,
as a _bottleneck region_ $P_{k}$, and refer to a _spontaneous_ transition past
$P_{k}$. (This is also called the _ghost of a saddle-node bifurcation_ in
strogatz1994nonlinear .) If all the connections corresponding to edges in the
graph $G$ are excitable, then the system contains an excitable network
attractor. If all connections are spontaneous, then we typically see a
periodic orbit, although we do not prove this.

##  3 Examples for the smooth activation function

###  3.1 Two vertex graph

For our first example we consider the connected graph with two vertices and a
single edge joining them, that is $a_{12}=1$, and $a_{ij}=0$ for
$(i,j)\neq(1,2)$. We use bifurcation analysis to show that the transition
between spontaneous and excitable dynamics is caused by a saddle-node
bifurcation, and find an approximation to the location of the saddle-node
bifurcation in parameter space.

The two-dimensional system of equations is:

$\displaystyle\dot{y}_{1}$
$\displaystyle=-y_{1}+w_{s}\phi(y_{1})+w_{m}\phi(y_{2}),$ (14)
$\displaystyle\dot{y}_{2}$
$\displaystyle=-y_{2}+w_{s}\phi(y_{2})+w_{p}\phi(y_{1}).$ (15)

In the $J_{i}$ variables, this becomes

$\displaystyle\dot{J}_{1}$
$\displaystyle=\frac{1}{\epsilon}J_{1}(1-J_{1})\left(-\phi^{-1}(J_{1})+w_{s}J_{1}+w_{m}J_{2}\right),$
(16) $\displaystyle\dot{J}_{2}$
$\displaystyle=\frac{1}{\epsilon}J_{2}(1-J_{2})\left(-\phi^{-1}(J_{2})+w_{s}J_{2}+w_{p}J_{1}\right),$

where $(J_{1},J_{2})\in(0,1)^{2}$. We note the following properties of the
function $g:(0,1)\rightarrow\mathbb{R}$, with $g(x)=\phi^{-1}(x)-w_{s}x$:

$g^{\prime}(x)=\frac{\epsilon}{x(1-x)}-w_{s},\quad
g^{\prime\prime}(x)=\frac{\epsilon(2x-1)}{x^{2}(1-x)^{2}},$
$\lim_{x\rightarrow 0}g(x)=-\infty,\quad\lim_{x\rightarrow 1}g(x)=\infty,$
$g\left(\frac{1}{2}\right)=\theta-\frac{w_{s}}{2},\quad
g^{\prime}\left(\frac{1}{2}\right)=4\epsilon-w_{s}.$

If $w_{s}>4\epsilon$, then $g$ has local extrema at $x_{+}$ and $x_{-}$, where

$x_{\pm}=\frac{1}{2}\pm\sqrt{\frac{1}{4}-\frac{\epsilon}{w_{s}}}.$ (17)

The $J_{1}$ and $J_{2}$ nullclines of system (16) are at, respectively

$J_{2}=\frac{1}{w_{m}}g(J_{1})\quad\text{and}\quad
J_{1}=\frac{1}{w_{p}}g(J_{2}).$ (18)

In figure 3, we show a sketch of the phase space of (16); the $J_{1}$
nullcline is shown in blue, and the $J_{2}$ nullcline in red. Solid dots show
stable equilibria, open dots show unstable equilibria and arrows show the
direction of flow. Note that we do not include any nullclines at $J_{i}=0$ or
$J_{i}=1$ because they are not in the domain of equations (16). As $w_{p}$ is
decreased (in the figures, moving from left to right), a saddle-node
bifurcation creates a pair of equilibrium solutions.

$J_{1}$$J_{1}$$J_{2}$$J_{2}$$x_{-}$$A$$B$ Figure 3:  Schematic illustration of
the phase space and nullclines of the system (16). The blue curves are the
$J_{1}$ nullclines, and the red curves are the $J_{2}$ nullclines. Solid dots
show stable equilibria and open dots show unstable equilibria. The parameter
$w_{p}$ is decreased from left figure to right figure, creating two further
equilibrium in a saddle-node bifurcation. The squares labelled $A$ and $B$ are
referred to in the proof of lemma 1.

###### Lemma 1

If $\theta<w_{s}$, and $0<\epsilon\ll 1$, then a saddle node bifurcation
occurs in system (16) when

$w_{p}=w_{p}^{SN}\equiv\epsilon\log\epsilon+\theta-\epsilon(1+\log
w_{s})+\frac{\epsilon^{2}}{w_{s}}+O(\epsilon^{3}).$ (19)

To begin the proof, we note that the saddle-node bifurcation will occur when
the points $A$ and $B$ (marked by squares in the left-hand panel of figure 3)
coincide. These points are defined as the intersection of the nullclines with
the line at $J_{2}=x_{-}=\frac{\epsilon}{w_{s}}+O(\epsilon^{2})$, i.e. at the
local extrema of the $J_{2}$ nullcline.

Let the $J_{1}$ coordinate of $A$ be

$\displaystyle J_{1}^{A}$ $\displaystyle=\frac{1}{w_{p}}g(x_{-})$
$\displaystyle=\frac{1}{w_{p}}\left(\epsilon\log\epsilon+\theta-\epsilon(1+\log
w_{s})+\frac{\epsilon^{2}}{w_{s}}\right)+O(\epsilon^{3})$

Let the $J_{1}$ coordinate of $B$ be $J_{1}^{B}$, and write
$J_{1}^{B}=1-\epsilon^{B}$, for some $\epsilon^{B}\ll 1$. Substituting this,
along with $J_{2}=x_{-}$, into the expression for the $J_{1}$ nullcline in
(18) gives

$x_{-}=\frac{1}{w_{m}}g(1-\epsilon^{B})$

Expanding in terms of the small quantities $\epsilon$ and $\epsilon^{B}$, this
gives

$\displaystyle\frac{\epsilon}{w_{s}}+O(\epsilon^{2})=$
$\displaystyle\frac{1}{w_{m}}\left(\theta-\epsilon(\log\epsilon^{B}+\epsilon^{B}+O({\epsilon^{B}}^{2}))-w_{s}+\epsilon^{B}w_{s}\right)$

which we rearrange to find

$\log\epsilon^{B}=\frac{\theta-
w_{s}}{\epsilon}+w_{s}\frac{\epsilon^{B}}{\epsilon}+O(1).$

Since we are assuming $\theta<w_{s}$, then $\epsilon^{B}$ is exponentially
small, that is $\epsilon^{B}=O(\epsilon^{n})$ for all $n\in\mathbb{N}$, thus,
$J_{1}^{B}=1+O(\epsilon^{n})$.

The points $A$ and $B$ collide when $J_{1}^{A}=J_{1}^{B}$, that is, when

$w_{p}=w_{p}^{SN}\equiv\theta+\epsilon\log\epsilon-\epsilon(1+\log
w_{s})+\frac{\epsilon^{2}}{w_{s}}+O(\epsilon^{3}).$

$\square$

For the default parameters (12), except for $w_{p}$, we find
$w_{p}^{SN}=0.3027$ (4 s.f.). Note that this means for $w_{p}=0.3$ we are
close to saddle node and there is an excitable connection with small
$\delta>0$. More generally, note that for any fixed $w_{s}$, as
$\epsilon\rightarrow 0$ we have $w_{p}^{SN}\rightarrow\theta$ as expected from
Theorem 2.1.

The following result gives an approximation of the positions of the equilibria
that are created in the saddle-node bifurcation. Methods similar to those used
in this proof are used in later sections for larger networks.

###### Lemma 2

If $\theta<w_{s}$, $0<\epsilon\ll 1$, and $0<\eta\ll\frac{\epsilon}{4}$, then
if $w_{p}=w_{p}^{SN}-\eta$, the system (16) has a pair of equilibria at

$(J_{1},J_{2})=\left(1,\frac{\epsilon}{w_{s}}\pm\frac{\sqrt{2\eta\epsilon}}{w_{s}}\right)+O(\epsilon^{2}).$

Recall that $x_{-}=\frac{\epsilon}{w_{s}}+O(\epsilon^{2})$. Thus,

$g^{\prime\prime}(x_{-})=-\frac{w_{s}^{2}}{\epsilon}\sqrt{1-\frac{4\epsilon}{w_{s}}}=-\frac{w_{s}^{2}}{\epsilon}+2w_{s}+O(\epsilon).$

Using the earlier results on the location of the $J_{1}$ nullcline, we will
have equilibria when

$\frac{g(J_{2})}{w_{p}}=1+O(\epsilon^{n}).$

Expanding $g$ about $J_{2}=x_{-}$ and writing $w_{p}=w_{p}^{SN}-\eta$ gives,

$\displaystyle g(J_{2})$
$\displaystyle=g(x_{-})+\frac{(J_{2}-x_{-})^{2}}{2}g^{\prime\prime}(x_{-})$
$\displaystyle~{}~{}~{}+O((J_{2}-x_{-})^{3}),$
$\displaystyle=w_{p}^{SN}-\eta+O(\epsilon^{n}),$
$\displaystyle(J_{2}-x_{-})^{2}$
$\displaystyle=-\frac{2\eta}{g^{\prime\prime}(x_{-})}+O(\epsilon^{3}),$

where the final line follows because $g(x_{-})=w_{p}^{SN}$. Substituting for
$g^{\prime\prime}(x_{-})$ then gives the result. $\square$

###  3.2 Three vertex cycle

Our second example is the cycle between three vertices shown schematically in
figure 4. As a heteroclinic cycle between equilibria, this system has been
studied extensively in the fields of populations dynamics ML75 , rotating
convection BH80  and symmetric bifurcation theory GH88 .

Figure 4:  Graph of the three vertex cycle.

We give some numerical examples of the dynamics of this system as realised by
the CTRNN excitable network, and use the continuation software AUTO auto  to
show that the transition from excitable to spontaneous dynamics occurs at a
saddle-node on an invariant circle (SNIC) bifurcation generating a periodic
orbit.

The deterministic equations realising this graph are:

$\displaystyle\dot{y}_{1}$
$\displaystyle=-y_{1}+w_{s}\phi(y_{1})+w_{m}\phi(y_{2})+w_{p}\phi(y_{3}),$
$\displaystyle\dot{y}_{2}$
$\displaystyle=-y_{2}+w_{s}\phi(y_{2})+w_{m}\phi(y_{3})+w_{p}\phi(y_{1}),$
(20) $\displaystyle\dot{y}_{3}$
$\displaystyle=-y_{3}+w_{s}\phi(y_{3})+w_{m}\phi(y_{1})+w_{p}\phi(y_{2}).$

We also consider the noisy case, using the setup given in equations (4).

Figure 5 show sample time series for two different parameter sets. On the
left, we show a noisy realisation with $w_{p}=0.3$, and on the right, a
periodic solution in the deterministic system (equations (1)) with
$w_{p}=0.305$.

$t$$t$$t$$t$$J_{k}$$J_{k}$$y_{k}$$y_{k}$ Figure 5: For the three-vertex cycle
(20), the top row of figures show time series of the $y_{k}$ variables in the
excitable (left) and spontaneous (right) cases. In the left column,
$w_{p}=0.3$ and $\sigma=0.05$. In the right column, $w_{p}=0.305$ and
$\sigma=0$. The bottom row show the $J_{k}$ variables. Other parameters are
$\epsilon=0.05$, $\theta=0.5$, $w_{s}=1$, $w_{m}=-0.7$.

Note that in both cases, for this system the $y_{k}$ variables oscillate
between three values: high ($y_{k}=Y_{A}=w_{s}=1$), intermediate
($y_{k}=Y_{L}=w_{p}=0.3$), and low ($y_{k}=Y_{T}=w_{m}=-0.7$), as the cells
shift between Active, Trailing and Leading. Only the first of these
corresponds to $J_{k}\approx 1$, as can be seen in the time series plots of
the $J_{k}$ variables in the lower panels of the figure, the other two
correspond to $J_{k}\approx 0$.

We compute a bifurcation analysis of the system (20) using the continuation
software AUTO auto97 .

SNIC$w_{p}$$T$$y_{1}$ Figure 6: For the three vertex cycle with equations
(20), the figure shows a bifurcation diagram as $w_{p}$ is varied. The top
panel shows the $y_{1}$ coordinate of equilibrium solutions, and the maximum
value of $y_{1}$ for periodic solutions. The lower panel shows the period of
the periodic solutions. Equilibrium solutions are shown by a thin line, and
periodic solutions by a thick line. Stable solutions are shown in red.
Bifurcation points are indicated by various shapes: Hopf bifurcations by
circles, saddle-node bifurcations by diamonds, saddle-node of periodic orbits
by triangles, and saddle-node on invariant circles (SNIC) by squares. Note
that there may be two squares for a single SNIC bifurcation because the
maximum value of $y_{1}$ on the periodic orbit is not the same as the value of
$y_{1}$ for the equilibria undergoing the saddle-node bifurcation. The SNIC
bifurcation labelled at $w_{p}=w_{p}^{SNIC}\approx 0.30287$ is the transition
between excitable and spontaneous dynamics.

Figure 6 shows a bifurcation diagram of this system as $w_{p}$ is varied.
Stable solutions are shown in red. There is a saddle-node on an invariant
circle (SNIC) bifurcation at $w_{p}=w_{p}^{SNIC}\approx 0.30287$. For
$w_{p}<w_{p}^{SNIC}$, the diagram shows a stable equilibrium solution with
$y_{1}\approx Y_{A}=1$ (and $y_{2}\approx Y_{L}$, $y_{3}\approx Y_{T}$). As
$w_{p}$ increases through $w_{p}^{SNIC}$, this equilibrium disappears in a
SNIC bifurcation creating a stable periodic orbit. Note that the period of the
periodic orbit asymptotes to $\infty$ as the SNIC bifurcation is approached.
Due to the symmetry, there are of course two further pairs of equilibria, one
pair with $y_{1}\approx Y_{T}$, $y_{2}\approx Y_{A}$, $y_{3}\approx Y_{L}$,
and another with $y_{1}\approx Y_{L}$, $y_{2}\approx Y_{T}$, $y_{3}\approx
Y_{A}$. The symmetry causes three saddle-node bifurcations to occur
simultaneously, creating the periodic orbit. If we were to instead choose the
$w_{p}$ to be different in each of the lines in (20), the saddle-nodes would
occur independently, and a periodic orbit would exist only if all three
$w_{p}$’s were greater than $w_{p}^{SNIC}$.

The time-series on the left-hand side of figure 5 has
$w_{p}=0.3<w_{p}^{SNIC}$. Without noise, at these parameter values, the system
would remain at one equilibrium point indefinitely. The noise acts as inputs
pushing the trajectory along the excitable connections. The time series on the
right-hand side of figure 5 has $w_{p}=0.305>w_{p}^{SNIC}$, and shows the
periodic orbit which has resulted from the SNIC bifurcation.

We note that this SNIC bifurcation occurs at approximately the same value of
$w_{p}$ as the saddle-node bifurcation found in section 3.1. This is not
surprising; using similar methods to those in the previous section, we can
show that to lowest order in $\epsilon$, the SNIC bifurcation occurs when
$w_{p}=w_{p}^{SN}$. For $w_{p}<w_{p}^{SNIC}$ there thus exists an excitable
network in the sense defined in appendix A.

###  3.3 Four node Kirk–Silber network

For our next example, we consider a graph with the structure of the
Kirk–Silber network KS94 , shown schematically in figure 7. This graph has one
vertex which has two outgoing edges, and the dynamics here are somewhat
different to vertices with only one outgoing edge. The bulk of this section is
devoted to an analysis of these differences, in particular, the possibility of
an additional equilibrium in the network attractor with two active cells.

$2$$1$$3$$4$ Figure 7:  Graph of the four vertex Kirk–Silber network.

The corresponding deterministic equations for this network are (moving
immediately into the $J_{i}$ variables):

$\displaystyle\dot{J}_{1}$
$\displaystyle=\frac{1}{\epsilon}J_{1}(1-J_{1})\left(-\phi^{-1}J_{1}+w_{s}J_{1}+w_{m}J_{2}+w_{p}J_{3}+w_{p}J_{4}\right),$
(21) $\displaystyle\dot{J}_{2}$
$\displaystyle=\frac{1}{\epsilon}J_{2}(1-J_{2})\left(-\phi^{-1}J_{2}+w_{p}J_{1}+w_{s}J_{2}+w_{m}J_{3}+w_{m}J_{4}\right),$
$\displaystyle\dot{J}_{3}$
$\displaystyle=\frac{1}{\epsilon}J_{3}(1-J_{3})\left(-\phi^{-1}J_{3}+w_{m}J_{1}+w_{p_{3}}J_{2}+w_{s}J_{3}+w_{t}J_{4}\right),$
$\displaystyle\dot{J}_{4}$
$\displaystyle=\frac{1}{\epsilon}J_{4}(1-J_{4})\left(-\phi^{-1}J_{4}+w_{m}J_{1}+w_{p_{4}}J_{2}+w_{t}J_{3}+w_{s}J_{4}\right).$

We can break the symmetry between $J_{3}$ and $J_{4}$ by choosing
$w_{p_{3}}\neq w_{p_{4}}$. In fact, in what follows, we will frequently write
$w_{p_{3}}=w_{p_{4}}+\Delta w$, for $\Delta w>0$, and choose $w_{p_{4}}=w_{p}$
for simplicity.

We consider first the dynamics near each of the vertices that have exactly one
outgoing edge (vertices 1, 3 and 4 in the graph; see figure 7). Again using
the same techniques that were used in section 3.1 (lemma 2), we can show that
for $w_{p},w_{m},w_{t}<\theta<w_{s}$, and $w_{p}=w_{p}^{SN}+\eta$,
$0<\eta<\frac{\epsilon}{4}$, there exist equilibria solutions at, for example

$(J_{1},J_{2},J_{3},J_{4})=\left(1,\frac{\epsilon}{w_{s}}\pm\frac{\sqrt{2\eta\epsilon}}{w_{s}},0,0\right)+O(\epsilon^{2}).$

That is, there is a transition from excitable to spontaneous dynamics (in this
case between cells $1$ and $2$, but also between cells $3$ and $1$ and cells
$4$ and $1$) as $w_{p}$ is increased through $w_{p}^{SN}$.

The dynamics close to the vertex with two outgoing connections (vertex 2) is
modified by the presence of the additional parameter $w_{t}$. Consider the
three dimensional subsystem of (21) with $J_{1}=0$, that is:

$\displaystyle\dot{J}_{2}$
$\displaystyle=\frac{1}{\epsilon}J_{2}(1-J_{2})\left(-\phi^{-1}J_{2}+w_{s}J_{2}+w_{m}J_{3}+w_{m}J_{4}\right)$
(22) $\displaystyle\dot{J}_{3}$
$\displaystyle=\frac{1}{\epsilon}J_{3}(1-J_{3})\left(-\phi^{-1}J_{3}+w_{p_{3}}J_{2}+w_{s}J_{3}+w_{t}J_{4}\right)$
$\displaystyle\dot{J}_{4}$
$\displaystyle=\frac{1}{\epsilon}J_{4}(1-J_{4})\left(-\phi^{-1}J_{4}+w_{p_{4}}J_{2}+w_{t}J_{3}+w_{s}J_{4}\right).$

We can perform a similar calculation to that shown in section 3.1 to show that
there is a section of the $J_{2}$ null-surface which lies asymptotically close
to the surface $J_{2}=1$. Equilibria solutions exist on this null-surface if
the $J_{3}$ and $J_{4}$ null-surfaces intersect there, that is, if there are
solutions to the pair of equations

$\displaystyle g(J_{3})$ $\displaystyle=w_{p_{3}}+w_{t}J_{4},$ (23)
$\displaystyle g(J_{4})$ $\displaystyle=w_{p_{4}}+w_{t}J_{3}.$ (24)

We assume without loss of generality that $w_{p_{3}}>w_{p_{4}}$ (i.e. $\Delta
w>0$) and then the arrangement of these curves is in one of the configurations
shown in figure 8. If $w_{t}<0$ (lower panel), equilibria solutions exist
(i.e. the red and blue curves intersect) for a range of $w_{p_{3}}$ and
$w_{p_{4}}$ with both larger than $w_{p}^{SN}$: that is, the transition to
spontaneous dynamics happens at a larger value of $w_{p_{j}}$ (than if
$w_{t}=0$). If $w_{t}>0$ (upper panel), the opposite happens: that is, the
transition to spontaneous dynamics occurs at a smaller value of $w_{p_{j}}$.
Solving these equations exactly requires solving a quartic equation, and the
resulting expression is not illuminating. We label the value of $w_{p}$ at
which this transition from spontaneous to excitable dynamics occurs as
$w_{p}^{SN^{\prime}}$, and note that this is a function of $w_{t}$, $w_{s}$,
$\epsilon$, $\theta$ as well as more generally, the number of Leading
directions from that cell.

$J_{3}$$J_{3}$$J_{4}$$J_{4}$ Figure 8:  Schematic illustration of the
nullclines of the system (22) system along the surface $J_{2}=1$. The blue
curves are the $J_{3}$ nullclines, and the red curves are the $J_{4}$
nullclines. In the upper panel, $w_{t}=0.05$, $w_{p_{3}}=0.30$,
$w_{p_{4}}=0.298$. In the lower panel, $w_{t}=-0.05$, $w_{p_{3}}=0.306$,
$w_{p_{4}}=0.304$. The dashed lines show the where the nullcline would lie if
$w_{t}=0$.

For the specific system (21), with $w_{p_{3}}=w_{p_{4}}+\Delta w=w_{p}+\Delta
w$, we thus have two conditions. If
$w_{p}<\min(w_{p}^{SN},w_{p}^{SN^{\prime}}-\Delta w)$ then the system is
excitable along all connections. If
$w_{p}>\max(w_{p}^{SN},w_{p}^{SN^{\prime}}-\Delta w)$ then a periodic orbit
will exist. If neither of these conditions holds then we will see excitable
connections in some places and spontaneous transitions in others.

(a)(b)(c)(d)$t$$t$$t$$t$$y_{j}$$y_{j}$$y_{j}$$y_{j}$ Figure 9: The figures
show time series of simulations of the $y_{j}$ variables in the Kirk–Silber
type network (21). In (a), $w_{p}=0.305$, and $\sigma=0$; in (b), $w_{p}=0.3$
and $\sigma=0.05$; in (c) $w_{p}=0.315$, and $\sigma=0$; in (d), $w_{p}=0.315$
and $\sigma=0.05$.

In figure 9 we show some example time-series of the system (21) (in the
$y_{k}$ coordinates). In panel (a), parameters are such that
$w_{p}>\max(w_{p}^{SN},w_{p}^{SN^{\prime}}-\Delta w)$, so we see a periodic
solution in the deterministic system. Note that the $y_{3}$ (yellow)
coordinate becomes close to $Y_{A}=1$ during this trajectory, but the $y_{4}$
(purple) coordinate does not: it switches between $Y_{L}=0.3$, $Y_{D}=0$ and
$Y_{T}=0.7$. In panel (b), parameters are such that
$w_{p}<\min(w_{p}^{SN},w_{p}^{SN^{\prime}}-\Delta w)$, so without noise the
trajectory would remain at a single equilibrium solution. Here we add noise
with $\sigma=0.05$, and the trajectory can be seen exploring the network. Note
that there are some transitions between $\xi_{2}$ ($y_{2}$ is red) and
$\xi_{3}$ ($y_{3}$ is yellow), and some from $\xi_{2}$ to $\xi_{4}$ ($y_{4}$
is purple). In panels (c) and (d), we increase $w_{p}$ further away from the
saddle-node bifurcation (further into the regime of spontaneous transitions)
and observe some qualitative differences in the trajectories. In the
deterministic case (c), the periodic solution now transitions from near the
bottleneck region $P_{2}$ to a region of phase space where $y_{3}$ and $y_{4}$
are both Active. We label this region of phase space as $P_{3,4}$. In the
noisy case (d) (which is also in the spontaneous regime), the trajectory makes
transitions from the bottleneck region $P_{2}$ to each of $P_{3}$, $P_{4}$,
and to $P_{3,4}$.

In the above simulations, we have used $w_{t}=0$, but we observe that the type
of qualitative behaviour observed depends both on the parameters $w_{p}$ and
$w_{t}$. Specifically, a sufficiently negative $w_{t}$ provides a suppression
effect, meaning that only a single cell $y_{k}$ can be active at any one time,
but the transitional value of $w_{t}$ depends on $w_{p}$. In figure 10 we show
maximum values of $y_{3}$ and $y_{4}$ along the periodic orbits as $w_{t}$ is
varied. It can be seen clearly here that the transition between periodic
orbits which visit $P_{3}$ (where $\max(y_{3})$ is significantly larger than
$\max(y_{4})$) and those which visit $P_{3,4}$ (where
$\max(y_{3})\approx\max(y_{4})$) is quite sharp.

$T$

$\max y_{3}$, $y_{4}$

$w_{t}$$w_{t}$ Figure 10: Behaviour of periodic orbits in the Kirk–Silber-type
four-node network (21) as the parameter $w_{t}$ is varied. The top panel shows
$\max(y_{3})$ (blue) and $\max(y_{4})$ (red) along the periodic orbit. The
bottom panel shows the period of the orbit on varying $w_{t}$.

We extend these results to show the behaviour as both parameters $w_{p}$ and
$w_{t}$ are varied in figure 11. The data in this figure shows the observed
behaviours for both noisy and deterministic systems as the parameters $w_{p}$
and $w_{t}$ are varied. The red lines are the curves $w_{p}=w_{p}^{SN}$
(dotted) and $w_{p}=w_{p}^{SN^{\prime}}$ (dashed). If $w_{p}$ is above both of
these lines then all transitions are spontaneous, and so a periodic orbit
exists in the system. If $w_{p}$ lies below either one (or both) of these
lines, then at least one of the transitions will be excitable and so there
will be no periodic solutions. The black line shows the boundary between those
periodic solutions which visit $P_{3}$ (to the left of the black line) and
those which visit $P_{3,4}$ (to the right of the black line), as determined by
the location of the sharp transition in calculations similar to those shown in
figure 10 for a range of $w_{p}$. The background colours are results from
noisy simulations. The colour indicates the ratio of transitions to $P_{3,4}$
to the total number of transitions to $P_{3}$, $P_{4}$ and $P_{3,4}$.
Interestingly, the noisy solutions require a much larger value of $w_{t}$ than
the deterministic ones to have a significant proportion of transitions to
$P_{3,4}$.

$w_{t}$$w_{p}$(b)(a)(c),(d) Figure 11:  Behaviour of the Kirk–Silber-type
four-node network (21) as the parameters $w_{t}$ and $w_{p}$ are varied. The
red lines are the curves $w_{p}=w_{p}^{SN}$ (dotted) and
$w_{p}=w_{p}^{SN^{\prime}}$ (dashed; determined numerically by solving a
quartic equation). For $w_{p}$ above both of these lines periodic solutions
exist in the deterministic system. To the left of the black line these
periodic solutions visit $P_{3}$, to the right they visit $P_{3,4}$. The
background colours are results from noisy simulations with $\sigma=0.05$. The
colour indicates the ratio of transitions to $P_{3,4}$ to the total number of
transitions to $P_{3}$, $P_{4}$ and $P_{3,4}$. The labelled dots give the
parameter values of the time-series plots in figure 9.

These changes in qualitative dynamics can be explained in terms of the three-
dimensional subsystem with $J_{1}=0$, given by equations (22). In this three-
dimensional system, there are stable equilibria at

$\displaystyle(J_{2},J_{3},J_{4})$ $\displaystyle=(0,1,0)+O(\epsilon)$
$\displaystyle(J_{2},J_{3},J_{4})$ $\displaystyle=(0,0,1)+O(\epsilon)$
$\displaystyle(J_{2},J_{3},J_{4})$ $\displaystyle=(0,1,1)+O(\epsilon)$

as well as further unstable/saddle equilibria. Recall that these equilibria
are not on the boundaries of the box (which are not part of the domain). In
figure 12 we show solutions from the full four-dimensional system (21)
projected onto the three-dimensional space with $J_{1}=0$. In panel (a), we
show the periodic solutions from the deterministic systems for five different
values of $w_{p}$, ranging from $w_{p}=0.309$ (left curve, dark purple), to
$w_{p}=0.3096$ (right curve, red) in increments of $0.0002$. It can be seen
that the first two of these trajectories approach the saddle equilibria
(marked as a blue dot) from one side of its stable manifold, and the latter
three from the other. The first three thus visit $P_{3}$, and the latter three
visit $P_{3,4}$. It is the transition of the periodic orbit across the stable
manifold which results in the rapid change in the qualitative behaviour of the
periodic orbit, and likely indicates that a homoclinic bifurcation to this
saddle point separates these behaviours. That is, the sharp transition in $T$
in figure 10 should actually extend to $\infty$ on both sides. In panel (b),
we show both noisy and deterministic trajectories with $w_{p}=0.315$. Note
that only one of the noisy trajectories follows the deterministic trajectory
closely: the majority visit $P_{3}$.

$J_{3}$$J_{2}$$J_{4}$$J_{3}$$J_{2}$$J_{4}$(a)(b) Figure 12:  The figures show
trajectories for the system (21) projected onto the three-dimensional space
with $J_{1}=0$. In panel (a), five periodic trajectories in the deterministic
system are shown for different values of $w_{p}$, ranging from $w_{p}=0.309$
(left curve, dark purple), to $w_{p}=0.3096$ (right curve, red) in increments
of $0.002$. In panel (b), we set $w_{p}=0.315$ (as in panels (c) and (d) of
figure 9). Noisy trajectories ($\sigma=0.05$) are shown in blue, and the
periodic orbit of the deterministic system is shown in purple. In both panels,
an equilibrium of the three-dimensional system (22) is shown by a blue dot.
Other parameters are $\varepsilon=0.05$, $\theta=0.5$, $w_{t}=0$, $w_{s}=1$,
$w_{m}=-0.5$. The arrows indicate the direction of flow.

###  3.4 A ten node network

In this section, we demonstrate the method of construction described in
section 2 for a larger network. Specifically, we randomly generated a directed
graph between 10 vertices, with the constraints that it contained no one-
loops, two-loops or $\Delta$-cliques, and such that the graph does not have
feed-forward structure (i.e. you cannot get ‘stuck’ in a subgraph by following
the arrows). The graph we consider is shown in figure 13.

Figure 13: A example of a directed graph between ten nodes with no loops of
order one or two and no $\Delta$-cliques. Three timeseries from realisations
of this as a network using CTRNNs are shown in figures 14-16.

We ran one simulation of the deterministic CTRNN system (1), and two
simulations of the noisy CTRNN system (4), and in each case, randomly
generated the entries for the $w_{p}$ in the equation (8). For the
deterministic system, the entries of $w_{p}^{ij}$ were chosen independently
from the uniform distribution $U(0.32,0.34)$. For the noisy systems, the
entries of $w_{p}^{ij}$ were chosen independently from the uniform
distribution $U(0.30,0.32)$. The remaining parameters were set at the default
parameter values given in (12) except $w_{t}=-0.3$ for the deterministic
system, and one of the noisy systems, and $\sigma=0.01$ for the noisy systems.
Note that for the deterministic parameter values there are bottlenecks in the
phase space close to the locations in phase space for the excitable states
that are present for the parameter values used in the stochastic cases.

$t$$y_{j}$ Figure 14: Trajectories of the CTRNN network for the directed graph
shown in figure 13, using the deterministic model (1). The top panel shows the
$y_{j}$ coordinates, where the colours correspond to the node colours in
figure 13. In the lower panel, each horizontal row corresponds to one node (as
labelled on the vertical axis), and the colour is blue when the corresponding
$J_{j}$ coordinate is close to zero, and yellow when it is close to one. That
is, the yellow segments indicate when each node is active. Parameters are as
described in the text.

The results of the simulations are shown in figures 14, 15 and 16. For the
deterministic simulation (figure 14), we see that the system has an attracting
period orbit, which visits the nodes in the order $1\rightarrow 4\rightarrow
7\rightarrow 3\rightarrow 8$. The entries of $w_{p}^{ij}$ were randomly
generated as described above, and we do not give them all here for space
reasons, but we note that in all cases in which a vertex in this cycle has two
‘choices’ for which direction to leave (in the graph shown in figure 13), the
attracting periodic orbit chooses the more unstable direction. That is, if $i$
is a vertex in the above cycle, and if $i\rightarrow j$ and $i\rightarrow k$
are connections in the directed graph, with $i\rightarrow j$ being a part of
the attracting periodic orbit, then $w_{p}^{ji}>w_{p}^{ki}$. Although we do
not prove here that this will always be the case, it is intuitively what one
might expect, that is, the connection from $i$ to $j$ is stronger than the
connection from $i$ to $k$.

$t$$y_{j}$ Figure 15: Trajectories of the CTRNN network for the directed graph
shown in figure 13, using the stochastic model (4). Lines and colours are
described in figure 14. Parameters are as described in the text, here
$w_{t}=-0.3$.

In the noisy simulation with $w_{t}=-0.3$ (figure 15), the equilibria are not
visited in a regular pattern, but random choices are made at each equilibria
from which there is more than one direction in which to leave. See, for
instance, the transition $3\rightarrow 8$ at $t\approx 75$, and the transition
$3\rightarrow 10$ at $t\approx 155$. The length of time spent near each
equilibria is also irregular; note for instance, the variable amount of time
spent near $\xi_{1}$ and $\xi_{3}$. In this simulation, because the transverse
parameter $w_{t}$ is sufficiently negative, only one node is active at any
given time.

By contrast, the simulation in figure 16 has $w_{t}=0$. Here, the transitions
again are made randomly, but without the suppression provided by the
transverse parameter it is possible to have multiple cells active at once. For
instance, at around $t=55$, both cells 1 and 5 become active at the same time;
they were both leading cells from the previously active cell 6. As cells 1
switches off, cell 4 becomes active, and as cell 5 switches off, cell 8
becomes active. The system continues to have two active cells around $t=250$,
at which point a third cell also becomes active. If the trajectory were to run
for longer then the number of active cells could decrease again, if an active
cell suppresses more than one previously active cells.

The entire excitable network attractor for this level of noise is clearly more
complicated than the design shown in figure 13, in that additional equilibria
(with more than one active cell) are accessible to those encoded and described
in Theorem 2.1. An interesting extension of this work would be to understand
which additional equilibria appear in a network attractor generated from a
given directed graph in this manner: figure 16 suggests that at least seven
levels of cell activity are needed to uniquely describe the states that can
appear when more than one cell becomes “active”. In particular, when a cell is
Trailing to more than one Active cell, the value of $y_{j}$ for that cell is
even lower than $Y_{T}$ (compare the top panel of figure 16 with the schematic
in figure 2).

$t$$y_{j}$ Figure 16: Trajectories of the CTRNN network for the directed graph
shown in figure 13, using the stochastic model (4). Lines and colours are
described in figure 14. Parameters are as described in the text, here
$w_{t}=0$.

##  4 Discussion

The main theoretical result of this paper is Theorem 2.1, which states that it
is possible to design the connection weight matrix of a CTRNN such that there
exists a network attractor with a specific graph topology embedded within the
phase space of the CTRNN. The graph topology is arbitrary except for minor
restrictions: namely there should be no loops of order one or two, and no
$\Delta$-cliques. Theorem 2.1 assumes a piecewise affine activation function,
but the examples in Section 3 suggest that the results generalise to CTRNN
using any suitable smooth activation function. More generally, note that the
coupled network is in some sense close to $N$ simultaneous saddle-node
bifurcations. However, the units are not weakly coupled and indeed this is
necessary to ensure that when one cell becomes active, the previous active
cell is turned off.

Theorem 2.1 proves the existence of an excitable network with threshold
$\delta$ where not only the connection weights, but also $\epsilon$ and
$\theta$ (properties of the activation function) may depend on $\delta$. We
believe that a stronger result will be true, namely that $\delta$ can be
chosen independent of properties of the activation function, and also that
this can be made an almost complete realisation by appropriate choice of
parameters.

###### Conjecture 1

Assume the hypotheses on the directed graph $G$ with $N$ vertices as in
Theorem 2.1 hold. Assume that $\epsilon>0$ is small and $\theta>0$. Then there
is a $\delta_{c}(\epsilon,\theta)>0$ such that for any $0<\delta<\delta_{c}$
there is an open set $\hat{W}_{\mathrm{ex}}\subset\mathbb{R}^{4}$. If the
parameters $(w_{s},w_{m},w_{t},w_{p})\in\hat{W}_{\mathrm{ex}}$ then the
dynamics of input-free equation (1) with $N$ cells and piecewise affine
activation function (3) and $w_{ij}$ defined by (8) contains an excitable
network attractor with threshold $\delta$ that gives an almost complete
realisation of the graph $G$.

The construction in Theorem 2.1 uses a comparatively sparse encoding of
network states - each of the $N$ vertices in the network is associated with
precisely one of the $N$ cells being in an active state. Indeed, the
connection weights (8) assign one of only four possible weights to each
connection, depending on whether that cell can become active next, was active
previously or neither. Other choices of weights will allow more dense
encoding: and many more than $N$ excitable states within a network of $N$
cells. However the combinatorial properties of the dynamics seem to be much
more difficult to determine and presumably additional connection weightings
will be needed, not the just four values considered in Theorem 2.1.

Section 3 illustrates specific examples of simple excitable networks for
smooth activation function (2) on varying parameters - this requires numerical
continuation to understand dependence on parameters even for fairly low
dimension. For this reason we expect that a proof of an analogous result to
Theorem 2.1 for the smooth activation function may be a lot harder. These
examples also give some insight into bifurcations that create the excitable
networks.

In general there is no reason that the realisation constructed in Theorem 2.1
is _almost complete_ (in the sense that almost all initial conditions in
$B_{\delta}(\xi_{k})$ evolve towards some $\xi_{l}$ with $a_{kl}=1$, analogous
to (ashwin2020almost, , Defn 2.6)). If it is not then other attractors may be
reachable from the excitable network. Conjecture 1 suggests that the
realisation can be made almost complete for small enough $\delta$: we expect
that for this we will require $w_{t}$ to be sufficiently negative. If $\delta$
is too large we cannot expect almost completeness: there are other stable
equilibria (notably the origin) that can be reached with large perturbations,
and the simulation in figure 16 shows that other equilibria may be reachable
from the network. It will be a challenge to strengthen our results to show
that the excitable network is an almost complete realization. However, the
examples studied in section 3 confirm that, at least for relatively simple
graphs, this conjecture is reasonable.

Our theoretical results are for networks with excitable connections. We expect
much of the behaviour described here is present in the spontaneous case, if
the coupling weights are chosen such that equilibria are replaced with
bottlenecks. In the absence of noise we expect to see a deterministic
switching between slow moving dynamics within bottlenecks. This dynamical
behaviour is very reminiscent of the stable heteroclinic channels described,
for example, in afraimovich2004heteroclinic ; rabinovich2020sequential .
However, stable heteroclinic switching models require structure in the form of
multiplicative coupling or symmetries that are not present in CTRNN or related
Wilson-Cowan neural models wilson1972excitatory  (see chow2019before  for a
recent review of related neural models). Other models showing sequential
excitation include chow2019before : this relies on a fast-slow decomposition
to understand various different modes of sequential activation in a neural
model of rhythm generation. It will be an interesting challenge to properly
describe possible output dynamics of our model in the case of bottlenecks.

We remark that asymmetry of connection weights is vital for constructing a
realisation as an excitable networks - indeed, the lack of two-cycles
precludes $a_{jk}=a_{kj}=1$. While this may be intuituively obvious, it was
not so obvious that we also need to exclude one-cycles and $\Delta$-cliques in
the graph to make robust realisations.

Finally, although we do not consider specific natural or machine learning
applications of CTRNN here, the structures found here may give insights that
give improved training for CTRNN. In particular, it seems plausible that CTRNN
may use excitable networks to achieve specific input-output tasks (especially
those requiring internal states). For example, recent work
ceni2020interpreting  demonstrates that echo state networks can create
excitable networks in their phase space to encode input-dependent behaviour.
It is also likely there are novel optimal training strategies that take
advantage of excitable networks, for example, choosing connection weights that
are distributed close to one of the four values we use.

###### Acknowledgements.

We gratefully acknowledge support from the Marsden Fund Council from New
Zealand Government funding, managed by Royal Society Te Apārangi. PA
acknowledges funding from EPSRC as part of the Centre for Predictive Modelling
in Healthcare grant EP/N014391/1. CMP is grateful for additional support from
the London Mathematical Laboratory.

##  Appendix A Definition of excitable network

This appendix (which extends ideas in Ashwin2016a ) gives formal definitions
for excitable networks considered in this paper. We say a system has an
_excitable connection for amplitude_ $\delta>0$ from one equilibrium $\xi_{i}$
to another $\xi_{j}$ if

$B_{\delta}(\xi_{i})\cap W^{s}(\xi_{j})\neq\emptyset,$

(where $B_{\delta}(\xi)$ is the open ball of radius $\delta$ centred at $\xi$)
and this connection has _threshold_ $\delta_{th}$ if

$\delta_{th}=\inf\\{\delta>0:B_{\delta}(\xi_{i})\cap
W^{s}(\xi_{j})\neq\emptyset\\}.$

We define the _excitable network of amplitude_ $\delta>0$ between the
equilibria $E=\\{\xi_{i}\\}$ to be the set

$\Sigma_{E}=\bigcup_{i,j=1}^{N}\\{\Phi_{t}(x):x\in
B_{\delta}(\xi_{i}),t>0\\}\cap W^{s}(\xi_{j}).$

We say the excitable network $\Sigma_{E}$ for amplitude $\delta$ _realises a
graph_ $G$ if each vertex $v_{i}$ in $G$ corresponds to an equilibrium
$\xi_{i}$ in $\Sigma_{E}$ and there is an edge in $G$ from $v_{i}$ to $v_{j}$
if only if there is a connection in $\Sigma_{E}$ for amplitude $\delta$ from
$\xi_{i}$ to $\xi_{j}$.

##  Appendix B Proof of Theorem 2.1

As in equations (9) and (10), for any $\delta<\frac{1}{2}$ we choose

$\epsilon=\dfrac{\delta}{8},~{}\theta=\frac{1}{2},~{}w_{s}=1,~{}w_{t}=0,$ (25)

and then $w_{p}$ and $w_{m}$ are given by

$w_{p}=\theta-\dfrac{\delta}{2},\quad
w_{m}=-(w_{s}-\theta)-\dfrac{\delta}{2}.$ (26)

We define equilibria $\xi_{k}$ as in Section 2.1, where we write the $j$th
component of the equilibrium as $[\xi_{k}]_{j}$:

$[\xi_{k}]_{j}=\begin{cases}Y_{A}&\mbox{ if }j=k,\\\ Y_{L}&\mbox{ if
}a_{kj}=1,\\\ Y_{T}&\mbox{ if }a_{jk}=1,\\\ Y_{D}&\mbox{ if }a_{kj}=0\mbox{
and }a_{jk}=0,\end{cases}$ (27)

where

$\displaystyle Y_{A}$ $\displaystyle:=w_{s}=1,\
Y_{L}:=w_{p}=\theta-\frac{\delta}{2},$ $\displaystyle Y_{T}$
$\displaystyle:=w_{m}=-(w_{s}-\theta)-\frac{\delta}{2},\ Y_{D}:=w_{t}=0,$

are the values of the Active, Leading, Trailing and Disconnected components,
respectively. Note that

$Y_{T}<Y_{D}<Y_{L}<\theta-\epsilon/2<\theta+\epsilon/2<Y_{A}.$ (28)

As mentioned before, the hypotheses of Theorem 2.1 imply that this labelling
is well defined and Figure 2 shows how a transition from cell 1 active to cell
2 active will occur in a general network. It is simple to check that (27) is
an equilbrium solution of (1) with activation function (3). Moreover,
$\xi_{k}$ is linearly stable with $n$ eigenvalues $-1$. We define
$[\mathcal{J}_{k}]_{j}:=\phi_{P}([\xi_{k}]_{j})$ then note that (28) implies
that

$[\mathcal{J}_{k}]_{j}:=\delta_{kj},$

in terms of the Kronecker $\delta_{kj}$.

$y_{k}$$y_{l}$$\mathcal{R}_{kl}$$\mathcal{R}_{l}$$\mathcal{S}_{kl}$$\mathcal{R}_{k}$$\xi_{k}$$\zeta_{k,l}$$\xi_{l}$$w_{s}$$w_{p}$$w_{s}$$w_{m}$$0$$\theta-4\epsilon$$\theta-2\epsilon$$\theta$$\theta+2\epsilon$$\theta+4\epsilon$
Figure 17: Illustration of a connection (shown in red) with amplitude
$\delta>0$ from $\xi_{k}$ to $\xi_{l}$, projected in the coordinates $y_{k}$
and $y_{l}$. The points $\xi_{k}$, $\xi_{l}$ and $0$ are linear sinks. The
point $\zeta_{k,l}$ is within $\delta$ of $\xi_{l}$ and limits to $\xi_{l}$ in
forwards time. The grey areas are regions of width $4\epsilon$ centred around
$y_{k,l}=\theta$ where $\phi_{P}(y_{k})$ and $\phi_{P}(y_{l})$ are non-
constant: these contain saddle and other equilibria (not shown).

We consider two cases. Case 1 is where $a_{kl}=1$ and we expect to see a
connection from $\xi_{k}$ to $\xi_{l}$. Case 2 is $a_{kl}=0$ and we don’t
expect a connection.

###  B.1 Case 1

Suppose $a_{kl}=1$ (recall that the lack of $1$-cycles means that $k\neq l$),
and then we define two regions of phase space, which we label
$\mathcal{R}_{l}$ and $\mathcal{R}_{kl}$, as follows:

$\mathcal{R}_{l}=\left\\{\mathbf{y}\ |\ y_{j}>\theta+\frac{\delta}{4},\ j=l;\
y_{j}<\theta-\frac{\delta}{4}\ \text{otherwise}\right\\},$
$\mathcal{R}_{kl}=\left\\{\mathbf{y}\ |\ y_{j}>\theta+\frac{\delta}{4},\
j=l,k;\ y_{j}<\theta-\frac{\delta}{4}\ \text{otherwise}\right\\}.$

The regions and the dynamics within them are shown schematically in figure 17.

Within $\mathcal{R}_{kl}$, if $a_{kl}=1$ then the dynamics are governed by the
equations

$\dot{y}_{j}=f_{j}^{(b)}:=-y_{j}+\phi^{(b)}_{j},$

where

$\phi_{j}^{(b)}:=\begin{cases}w_{s}+w_{m}&\mbox{ if }j=k,\text{(AT)}\\\
w_{s}+w_{p}&\mbox{ if }j=l,\text{(LA)}\\\ w_{p}&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(1,0,0,0),\text{(LD)}\\\ w_{m}&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(0,1,0,0),\text{(TD)}\\\ w_{p}&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(0,0,1,0),\text{(DL)}\\\ w_{m}&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(0,0,0,1),\text{(DT)}\\\ w_{m}+w_{p}&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(0,1,1,0),\text{(TL)}\\\ 0&\mbox{ if
}(a_{kj},a_{jk},a_{lj},a_{jl})=(0,0,0,0),\text{(DD)}\end{cases}$ (29)

where the type of coordinate (see figure 2) is given in parentheses on each
line. Within $\mathcal{R}_{l}$, the dynamics are governed by the equations

$\dot{y}_{j}=f_{j}^{(c)}:=-y_{j}+\phi^{(c)}_{j},$ (30)

where

$\phi_{j}^{(c)}:=\begin{cases}w_{m}&\mbox{ if }j=k,\text{(AT)}\\\ w_{s}&\mbox{
if }j=l,\text{(LA)}\\\ w_{p}&\mbox{ if }a_{lj}=1\mbox{ and }j\neq
l,\text{(DL/TL)}\\\ w_{m}&\mbox{ if }a_{jl}=1,\text{(DT)}\\\ 0&\mbox{ if
}a_{lj}=0\mbox{ and }a_{jl}=0,\text{(LD/TD/DD)}\end{cases}$ (31)

The equilbrium $\xi_{l}$ lies in the interior of the region $\mathcal{R}_{l}$,
whereas $\xi_{k}$ lies in the interior of the region $\mathcal{R}_{k}$. We
show that there is a connection of amplitude $\delta$ from $\xi_{k}$ to any
$\xi_{l}$ with $a_{kl}=1$, by considering a trajectory starting at

$\zeta_{k,l}=\xi_{k}+\delta e_{l},$

where $e_{l}$ is a unit vector in the $l$-direction: clearly
$|\zeta_{k,l}-\xi_{k}|=\delta$ (see figure 17). Note that

$[\zeta_{k,l}]_{j}=\begin{cases}Y_{A}&\mbox{ if }j=k,\\\ Y_{L}+\delta&\mbox{
if }j=l,\\\ Y_{L}&\mbox{ if }a_{kj}=1\mbox{ and }j\neq l,\\\ Y_{T}&\mbox{ if
}a_{jk}=1,\\\ Y_{D}&\mbox{ if }a_{kj}=0\mbox{ and }a_{jk}=0,\end{cases}$ (32)

where $Y_{A}$ etc are given in equation (11), and
$Y_{L}+\delta=\theta+\delta/2>\theta$. We define
$[\mathcal{J}_{k,l}]_{j}:=\phi_{P}([\zeta_{k,l}]_{j})$, and then

$[\mathcal{J}_{k,l}]_{j}=\delta_{jk}+\delta_{jl}.$

Now, we note that $\zeta_{k,l}\in\mathcal{R}_{kl}$. We next show that a
trajectory with initial condition at $\zeta_{k,l}$ will asymptotically
approach $\xi_{l}$. We show first that the trajectory enters $\mathcal{R}_{l}$
in finite time. Then, since in $\mathcal{R}_{l}$, the flow is linear, with
stable equilibrium $\xi_{l}$, all trajectories in $\mathcal{R}_{l}$ eventually
approach $\xi_{l}$.

Define $\mathcal{S}_{kl}$ to be the region between $\mathcal{R}_{l}$ and
$\mathcal{R}_{kl}$, namely,

$\mathcal{S}_{kl}=\left\\{\mathbf{y}\ |\ \theta-\frac{\delta}{4}\leq
y_{k}\leq\theta+\frac{\delta}{4};y_{l}>\theta+\frac{\delta}{4};\
y_{j}<\theta-\frac{\delta}{4}\ \text{ for }j\neq k,l\right\\}.$

First consider the dynamics of all $y_{j}$, with $j\neq k,l$. This means that
$y_{j}<\theta-\delta/4$ for all points in
$\mathcal{T}_{kl}\equiv\mathcal{R}_{kl}\cup\mathcal{S}_{kl}\cup\mathcal{R}_{l}$.
While the trajectory remains in $\mathcal{T}_{kl}$, it can be shown that
$\dot{y}_{j}$ is negative along the line with $y_{j}=\theta-\delta/4$. Hence,
none of the $y_{j}$ will leave $\mathcal{T}_{kl}$.

Within $\mathcal{S}_{kl}$, the dynamics of $y_{k}$ and $y_{l}$ are governed by

$\displaystyle\dot{y}_{k}=$ $\displaystyle-y_{k}+w_{m}+\phi_{P}(y_{k}),$ (33)
$\displaystyle\dot{y}_{l}=$ $\displaystyle-y_{l}+w_{s}+w_{p}\phi_{P}(y_{k}).$
(34)

Consider the equation for $\dot{y}_{k}$ in $\mathcal{R}_{kl}$,
$\mathcal{S}_{kl}$ and $\mathcal{R}_{l}$, namely:

$\dot{y}_{k}=\begin{cases}-y_{k}+w_{s}+w_{m}&\text{if}\
\mathbf{y}\in\mathcal{R}_{kl},\\\ -y_{k}+\phi_{P}(y_{k})+w_{m}&\text{if}\
\mathbf{y}\in\mathcal{S}_{kl},\\\ -y_{k}+w_{m}&\text{if}\
\mathbf{y}\in\mathcal{R}_{l}.\end{cases}$

Recall that $w_{m}=-(w_{s}-\theta)-\frac{\delta}{2}$, and
$0\leq\phi_{P}(y_{k})\leq 1$. We can use these bounds to show that

$\dot{y}_{k}\leq\begin{cases}-\dfrac{3\delta}{4}&\text{if}\
\mathbf{y}\in\mathcal{R}_{kl},\\\ -\dfrac{\delta}{4}&\text{if}\
\mathbf{y}\in\mathcal{S}_{kl},\end{cases}$

Furthermore, if $\mathbf{y}\in\mathcal{R}_{l}$ and $y_{k}>0$, then
$\dot{y}_{k}<-w_{m}<0$. In particular, we note that for
$\mathbf{y}\in\mathcal{T}_{kl}$, with $y_{k}>0$, $\dot{y}_{k}$ is negative and
bounded below zero.

Now, consider the equation for $\dot{y}_{l}$ in $\mathcal{R}_{kl}$,
$\mathcal{S}_{kl}$ and $\mathcal{R}_{l}$, namely:

$\dot{y}_{l}=\begin{cases}-y_{l}+w_{s}+w_{p}&\text{if}\
\mathbf{y}\in\mathcal{R}_{kl},\\\ -y_{l}+w_{s}+w_{p}\phi_{P}(y_{k})&\text{if}\
\mathbf{y}\in\mathcal{S}_{kl},\\\ -y_{l}+w_{s}&\text{if}\
\mathbf{y}\in\mathcal{R}_{l},\end{cases}$

We use this to compute $\dot{y}_{l}$ along the lower boundary of the three
regions, $\mathcal{R}_{kl}$, $\mathcal{S}_{kl}$ and $\mathcal{R}_{l}$, that
is, the line $y_{l}=\theta+\frac{\delta}{4}$, and we find

$\dot{y}_{l}=\begin{cases}w_{s}-\frac{3\delta}{4}&\text{if}\
\mathbf{y}\in\mathcal{R}_{kl},\\\
w_{s}-\theta-\frac{\delta}{4}+\left(\theta-\frac{\delta}{2}\right)\phi_{P}(y_{k})>w_{s}-\frac{3\delta}{4}&\text{if}\
\mathbf{y}\in\mathcal{S}_{kl},\\\ w_{s}-\theta-\frac{\delta}{4}&\text{if}\
\mathbf{y}\in\mathcal{R}_{l}.\end{cases}$

Since $w_{s}>\frac{3\delta}{4}$ and $w_{s}>\theta+\frac{\delta}{4}$, we see
that $\dot{y}_{l}>0$ in all three cases.

Combining our knowledge of $\dot{y}_{l}$ and $\dot{y}_{k}$ tells us that a
trajectory which starts in $\mathcal{R}_{kl}$, or more specifically, a
trajectory starting in a small neighbourhood of $\zeta_{k,l}$ will have
monotonic decreasing $y_{k}$ component until (at least) $y_{k}=0$.
Furthermore, the $y_{l}$ component cannot decrease below
$y_{l}=\frac{1}{2}+\frac{\delta}{4}$. Thus the trajectory will move through
$\mathcal{S}_{kl}$ and into $\mathcal{R}_{l}$ in a bounded time.

Within $\mathcal{R}_{l}$, $\xi_{l}$ is a linearly stable fixed point. In
summary, we have demonstrated that if $a_{kl}=1$ then there is a connection
from a $\delta$-neighbourhood of $\xi_{k}$ to $\xi_{l}$. Moreover, as the
equilibria are linearly stable and having a connection is an open condition,
the realisation will persist for an open set of parameters.

###  B.2 Absence of excitable connections for edges absent from $G$

Now suppose that $a_{lk}=a_{kl}=0$. Then the dynamics for $x_{k}$ and $x_{l}$
is shown schematically in figure 18. Equilibria are shown with dots, and all
equilibria shown in this figure are linearly stable. Note that the equilibrium
$\xi_{k}$ has $y_{k}=Y_{A}=1$, and $y_{l}=Y_{D}=0$. It is clear that there are
no small perturbations which allow for a connection between $\xi_{k}$ and
$\xi_{l}$.

For $\theta=1/2$, if $\delta<1/4$ then for any $k$ and $j\neq k$ such that
$a_{kj}=0$, all trajectories starting in

$\zeta_{k,l}=\xi_{k}+ae_{k}+be_{l},$

with $a^{2}+b^{2}\leq\delta^{2}$. In particular, perturbations of the form
(32) will return to $\xi_{k}$. This is because this set is remains in the
region of validity for the equivalent of (30) for which the only attractor is
$\xi_{k}$. In the case $a_{kl}=0$ and $a_{lk}=1$ a similar argument holds as
the phase portrait corresponds to figure 17 reflected in the diagonal.

This shows that no perturbations of amplitude $\delta$ within the
$(x_{k},x_{l})$ plane that give a connection from $\xi_{k}$ to $\xi_{l}$ of
amplitude $\delta$. However, we cannot rule out the existence of connections
from other locations in $B_{\delta}(\xi_{k})$ to $\xi_{l}$. This would be
needed to prove that the realisation is almost complete.

$y_{k}$$y_{l}$$\xi_{k}$$\xi_{l}$$w_{s}$$w_{s}$$0$$\theta-2\epsilon$$\theta$$\theta+2\epsilon$
Figure 18: Schematic diagram showing the dynamics in the $y_{k}$-$y_{l}$ plane
when $a_{kl}=a_{lk}=0$. Equilibria are shown with dots, and all equilibria are
linearly stable.

## References

  * (1) Afraimovich, V., Zhigulin, V., Rabinovich, M.: On the origin of reproducible sequential activity in neural circuits.  Chaos: An Interdisciplinary Journal of Nonlinear Science 14(4), 1123–1129 (2004)
  * (2) Afraimovich, V.S., Rabinovich, M.I., Varona, P.: Heteroclinic contours in neural ensembles and the winnerless competition principle.  International Journal of Bifurcation and Chaos 14(04), 1195–1208 (2004)
  * (3) Ashwin, P., Castro, S.B., Lohse, A.: Almost complete and equable heteroclinic networks.  Journal of Nonlinear Science 30(1), 1–22 (2020)
  * (4) Ashwin, P., Postlethwaite, C.: Designing Heteroclinic and Excitable Networks in Phase Space Using Two Populations of Coupled Cells.  J. Nonlinear Science 26(2), 345–364 (2016).  DOI 10.1007/s00332-015-9277-2
  * (5) Ashwin, P., Postlethwaite, C.: Sensitive finite-state computations using a distributed network with a noisy network attractor.  IEEE Transactions on Neural Networks and Learning Systems 29(12), 5847–5858 (2018)
  * (6) Beer, R.D.: On the dynamics of small continuous-time recurrent neural networks.  Adaptive Behavior 3(4), 469–509 (1995)
  * (7) Bhowmik, D., Nikiforou, K., Shanahan, M., Maniadakis, M., Trahanias, P.: A reservoir computing model of episodic memory.  In: 2016 International Joint Conference on Neural Networks (IJCNN), pp. 5202–5209. IEEE (2016)
  * (8) Blynel, J., Floreano, D.: Exploring the T-maze: Evolving learning-like robot behaviors using CTRNNs.  In: Applications of Evolutionary Computing, pp. 593–604. Springer Berlin Heidelberg, Berlin, Heidelberg (2003)
  * (9) Busse, F.M., Heikes, K.E.: Convection in a rotating layer: A simple case of turbulence.  Science 208, 173–175 (1980)
  * (10) Ceni, A., Ashwin, P., Livi, L.: Interpreting recurrent neural networks behaviour via excitable network attractors.  Cognitive Computation 12(2), 330–356 (2020)
  * (11) Chow, C.C., Karimipanah, Y.: Before and beyond the Wilson–Cowan equations.  Journal of Neurophysiology 123(5), 1645–1656 (2020).  DOI 10.1152/jn.00404.2019
  * (12) Doedel, E.J., Champneys, A.R., Dercole, F., Fairgrieve, T.F., Kuznetsov, Y.A., Oldeman, B., Paffenroth, R., Sandstede, B., Wang, X., Zhang, C.: AUTO-07P: Continuation and bifurcation software for ordinary differential equations (2007)
  * (13) Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Sandstede, B., Wang, X.J.: AUTO97: Continuation and bifurcation software for ordinary differential equations.  Tech. rep., Department of Computer Science, Concordia University, Montreal, Canada (1997).  (Available by FTP from ftp.cs.concordia.ca in directory pub/doedel/auto)
  * (14) Funahashi, K., Nakamura, Y.: Approximation of dynamical systems by continuous time recurrent neural networks.  Neural Networks 6, 801–806 (1993)
  * (15) Gouzé, J.L., Sari, T.: A class of piecewise linear differential equations arising in biological models.  Dynamical Systems 17, 299–316 (2003)
  * (16) Guckenheimer, J., Holmes, P.: Structurally stable heteroclinic cycles.  Math. Proc. Camb. Phil. Soc. 103, 189–192 (1988)
  * (17) Harris, J., Ermentrout, B.: Bifurcations in the Wilson–Cowan equations with nonsmooth firing rate.  SIAM Journal on Applied Dynamical Systems 14(1), 43–72 (2015)
  * (18) Hopfield, J.J., Tank, D.W.: Neural computation of decisions in optimization problems.  Biological Cybernetics 52(3), 141–152 (1985)
  * (19) Hutt, A., beim Graben, P.: Sequences by metastable attractors: Interweaving dynamical systems and experimental data.  Frontiers in Applied Mathematics and Statistics 3, 11 (2017).  DOI 10.3389/fams.2017.00011
  * (20) Kirk, V., Silber, M.: A competition between heteroclinic cycles.  Nonlinearity 7, 1605–1621 (1994)
  * (21) Manjunath, G., Tiño, P., Jaeger, H.: Theory of input driven dynamical systems.  ESANN 2012 proceedings (2012)
  * (22) May, R., Leonard, W.: Nonlinear aspects of competition between three species.  SIAM J. Appl. Math. 29, 243–253 (1975)
  * (23) Nikiforou, K.: The dynamics of continuous-time recurrent neural networks and their relevance to episodic memory.  Ph.D. thesis, Imperial College, London (2019)
  * (24) Rabinovich, M., Volkovskii, A., Lecanda, P., Huerta, R., Abarbanel, H., Laurent, G.: Dynamical encoding by networks of competing neuron groups: winnerless competition.  Physical review letters 87(6), 068102 (2001)
  * (25) Rabinovich, M.I., Huerta, R., Varona, P., Afraimovich, V.S.: Generation and reshaping of sequences in neural systems.  Biological cybernetics 95(6), 519–536 (2006)
  * (26) Rabinovich, M.I., Zaks, M.A., Varona, P.: Sequential dynamics of complex networks in mind: Consciousness and creativity.  Physics Reports (2020)
  * (27) Strogatz, S.H.: Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering.  Addison-Wesley (1994)
  * (28) Tuci, E., Quinn, M., Harvey, I.: An evolutionary ecological approach to the study of learning behavior using a robot-based model.  Adaptive Behavior 10(3-4), 201–221 (2002)
  * (29) Wilson, H.R., Cowan, J.D.: Excitatory and inhibitory interactions in localized populations of model neurons.  Biophysical Journal 12(1), 1–24 (1972)
  * (30) Yamauchi, B.M., Beer, R.D.: Sequential behavior and learning in evolved dynamical neural networks.  Adaptive Behavior 2(3), 219–246 (1994)