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

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

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.

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$.

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

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

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

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.

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:

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

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$:

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

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

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.

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

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

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

which we rearrange to find

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

$\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.

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

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

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

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

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 .

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:

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$.

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 .

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.

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.

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

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

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:

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

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.

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.

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.

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}$.

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