Royal Society Publishing

Delays in activity-based neural networks

Stephen Coombes, Carlo Laing


In this paper, we study the effect of two distinct discrete delays on the dynamics of a Wilson–Cowan neural network. This activity-based model describes the dynamics of synaptically interacting excitatory and inhibitory neuronal populations. We discuss the interpretation of the delays in the language of neurobiology and show how they can contribute to the generation of network rhythms. First, we focus on the use of linear stability theory to show how to destabilize a fixed point, leading to the onset of oscillatory behaviour. Next, we show for the choice of a Heaviside nonlinearity for the firing rate that such emergent oscillations can be either synchronous or anti-synchronous, depending on whether inhibition or excitation dominates the network architecture. To probe the behaviour of smooth (sigmoidal) nonlinear firing rates, we use a mixture of numerical bifurcation analysis and direct simulations, and uncover parameter windows that support chaotic behaviour. Finally, we comment on the role of delays in the generation of bursting oscillations, and discuss natural extensions of the work in this paper.


1. Introduction

Delays arise naturally in models of neurobiological systems. For example, the finite speed of an action potential (AP) propagating along an axon means that spike signalling between neurons depends upon how far apart they are. Hence, the interest in understanding network models with space-dependent delays, as in Laing & Coombes (2006). Upon arrival at a synaptic contact point, the transduction of an electrical signal into a biochemical signal and back again, to a post-synaptic potential (PSP), gives rise to a further delay. Yet another delay is associated with the spread of the PSP through the dendritic tree of the neuron to the cell body, where further APs can be initiated. It is now quite common to model both these forms of signal processing using either a form of distributed delay, as in Laing & Longtin (2003), or as a simple fixed or discrete delay. For an excellent review of the role of time delays in neural systems, we refer the reader to the paper by Campbell (2007). The effects of such delays can be quite varied. Although commonly associated with the generation of oscillations (Plant 1981), delays can also lead to oscillator death (Reddy et al. 1998), control phase locking (Coombes & Lord 1997) and underlie multi-stability (Shayer & Campbell 2000).

In this paper, we focus on the dynamics of two-population neural models with the incorporation of two discrete delays. In particular, we will work with the well-known Wilson & Cowan (1972) model. Such activity-based models are expected to provide a caricature of the behaviour of more realistic spiking networks when the time scale of synaptic processing is much longer than the membrane time constant of a typical cell (Ermentrout 1998). This is perhaps most clearly demonstrated by recent work of Roxin et al. (2005), which further emphasizes that a single delay in the activity-based representation can further improve the match with spiking networks. The delay in the activity-based model is interpreted by them as describing the time course of the AP initiation. However, an alternative interpretation of this delay is that it is necessary to adequately model the time lag involved in generating a rate-based representation of a spiking network. In particular, single-neuron firing rates (for slow synapses) will be largely determined by the steady-state values of non-spiking currents, and thus the delay may be more naturally interpreted as the time for these currents to relax. In any case, this paper will show how to analyse a delayed neural network with a hybrid approach, combining linear stability theory, the construction of periodic orbits (for piecewise constant nonlinear firing rate functions) and numerical techniques.

2. The model

As discussed earlier, under certain approximations, spiking network models can be reduced to just a few variables. One famous example is the Wilson & Cowan (1972) model that describes the evolution of a network of synaptically interacting neuronal populations (typically one being excitatory and the other inhibitory). In the presence of delays, this model takes the formEmbedded Image(2.1)Here, u and v represent the synaptic activity of the two populations, with a relative time scale for the response set by α−1. The architecture of the network is fixed by the weights a, b, c, d, while θu,v describe background drives (biases). The firing rate function f is commonly chosen as a sigmoidEmbedded Image(2.2)which satisfies the (Riccati) equation f′=βf(1−f), with β>0. The fixed delays τ1 and τ2 distinguish between the delayed self-interactions and the delayed cross-interactions. The delay differential equation (DDE) model (2.1) is similar, although not equivalent, to voltage-based models, which have linear combinations of sigmoidal functions of the different variables on the right-hand side (Marcus & Westervelt 1989; Olien & Bélair 1997; Wei & Ruan 1999; Shayer & Campbell 2000). Restrictions of the parameter choices recover a number of models already considered in the literature, such as that of (i) Glass et al. (1988), when a<0, b=0, (ii) Chen & Wu (1999), when α=1, a=d=0 and b=c>0, and (iii) Battaglia et al. (2007), when α=1, a=d<0, b=c>0 and f(z) is a threshold-linear firing rate.

Before analysing the full DDE system, it is first useful to describe the dynamics in the absence of delays, where we recover the basic Wilson and Cowan model. For τ1=τ2=0, it is straightforward to find the values of θu and θv corresponding to Hopf and saddle-node bifurcations. The point (u*, v*) is an equilibrium if there is a solution to the pair of equationsEmbedded Image(2.3)where Embedded Image. The Jacobian matrix is thereforeEmbedded Image(2.4)Thus, the conditions for a Hopf bifurcation (HB) areEmbedded Image(2.5)Eliminating v* asEmbedded Image(2.6)we can then plot the fixed point equation (parametrically) in the (θu, θv) plane, as in figure 1. A similar procedure can be used to determine the locus of saddle-node (SN) bifurcations defined by det L=0, as well as the Bogdanov–Takens bifurcation defined by det L=0 and Tr L=0 (when the SN and HB curves intersect). Indeed the Wilson and Cowan model also supports a saddle node on an invariant circle bifurcation (when the SN curve lies between the two HB curves), and can also support a saddle–separatrix loop and a double limit cycle. See Hoppensteadt & Izhikevich (1997, ch. 2) for a detailed discussion.

Figure 1

Hopf (HB, dashed line) and saddle-node (SN, solid line) bifurcation sets in the Wilson–Cowan network (no delays) with a mixture of excitatory and inhibitory connections for α=1, a=−b=c=10, d=2 and β=1.

3. Linear stability analysis of fixed point

The existence of an equilibrium is, of course, independent of any delays. Many authors have described in detail how the presence of delays affects the stability of an equilibrium, and here we follow the spirit of work by Marcus & Westervelt (1989), Wei & Ruan (1999) and Giannakapoulos & Zapp (2001). In the presence of delays, the linearized equations of motion have solutions of the form Embedded Image. Demanding that the amplitudes Embedded Image be non-trivial gives a condition on λ that may be written in the form Embedded Image(λ)=0, whereEmbedded Image(3.1)and the equilibrium (u*, v*) is given by the simultaneous solution of (2.3). For λEmbedded Image, we see that λ=0 whenEmbedded Image(3.2)where κ1=aβu*(1−u*), κ2=dβu*(1−v*) and κ3=bcβ2u*v*(1−u*)(1−v*). Thus, a real instability of a fixed point is defined by (3.2) and is independent of (τ1, τ2). Referring back to the analysis of §2, we see that this is identical to the condition for a saddle-node bifurcation. By contrast, a dynamic instability will occur whenever λ=iω for ω≠0, where ωEmbedded Image. The bifurcation condition in this case is defined by the simultaneous solution of the equations Re Embedded Image(iω)=0 and Im Embedded Image(iω)=0, namelyEmbedded Image(3.3)Embedded Image(3.4)For the parameters that ensure ω≠0, we shall say that the simultaneous solution of equations (3.3) and (3.4) defines a Hopf bifurcation at Embedded Image. More correctly, we should also ensure that as the delays pass through this critical point that the rate of change of Re λ is non-zero (transversality) and that there are no other eigenvalues with zero real part (non-degeneracy).

Interestingly, the models with two delays can lead to an interference effect whereby although either delay, if long enough, can bring about instability, there is a window of (τ1, τ2) where the solutions are stable to Hopf bifurcations. This is nicely discussed in the book by MacDonald (1989, ch. 6); see also the book by Stépán (1989, ch. 3). An example of this effect, obtained by computing the locus of Hopf bifurcations according to the above prescription, is shown in figure 2. A similar figure, showing a band of stability that lies between two broad regions of instability, is found in the work of Murdoch et al. (1987).

Figure 2

A bifurcation diagram showing the stability of the equilibrium in the Wilson and Cowan model with two delays. The parameters are the same as shown in figure 1 with (θu, θv)=(−2, −4).

4. Synchronous and anti-synchronous solutions

In general, despite linear stability analysis showing where to look, it is a challenge to find periodic solutions in closed form. Moreover, determining their stability is a problem that, in general, is best examined with numerical tools. However, some results are known about the phase relationship between the two populations during an oscillation. In particular, Chen et al. (2000) have shown that for α=1, θu=θv, a=d=0 and b=c that every non-constant solution of (2.1) is either synchronous or phase locked. Here, we explore the explicit construction of such solutions in the limit of high gain (β→∞), so that f(z)=H(z), where H is the Heaviside step function. Such equations are commonly encountered in physiological control systems (Glass et al. 1988; Longtin & Milton 1998). For example, in figure 3, we show a coexisting synchronous and anti-synchronous stable periodic orbit in a network with purely inhibitory connections. Previous work on the analysis of periodic orbits in delayed neural networks with Heaviside nonlinearity can be found in Guo et al. (2005).

Figure 3

Coexisting (a) synchronous and (b) anti-synchronous solutions for f(z)=H(z). The parameters are α=1, a=d=−1, b=c=−0.4, θu=θv=0.7, τ1=1 and τ2=1.4.

(a) Inhibitory network

We first consider a purely inhibitory network with a, b, c, d<0, with some bias θu=θv and α=1. Regarding a synchronous T-periodic solution, u(t)=v(t) with u(t+T)=u(t), like that shown in figure 3a, we parametrize such a solution in terms of two fundamental times T1,2 and the maxima and minima A± of the orbit. Here, T1 denotes the time spent on the decreasing part of the trajectory, and T2 that spent on the rising phase. Exploiting the piecewise linear nature of the dynamics, we then have thatEmbedded Image(4.1)Embedded Image(4.2)Embedded Image(4.3)Embedded Image(4.4)Solving these we obtain the period of oscillation T=T1+T2, whereEmbedded Image(4.5)and Embedded Image. The amplitude of the oscillation is A=A+A=(a+b+s)/s.

Similarly, to analyse an anti-synchronous solution, u(t)=v(t+T/2) with u(t)=u(t+T), as in figure 3b, we note that by symmetry, the rising and falling phases have the same duration, say T1. For the parameters considered, we find that τ1<T1<τ2, and we obtain the relationshipsEmbedded Image(4.6)Embedded Image(4.7)Embedded Image(4.8)Solving the above we find that T1 satisfies the transcendental equationEmbedded Image(4.9)The period T is 2T1 and the absolute amplitude of oscillation, A=A+A, is given byEmbedded Image(4.10)

(b) Excitatory self-feedback

For a single population with self-feedback, it is also possible to construct periodic solutions (for a Heaviside firing rate). Here, we consider just the evolution of u with a=1, b=0, θu=−h, h>0 and τ1=τd, a fixed delay. An example of a periodic trajectory is shown in figure 4. It is natural to parametrize the solution in terms of the four unknowns A± and T±, which denote the largest (A+) and smallest (A) values of the trajectory and the times spent above (T+) and below (T) the threshold h. The trajectory increases from A for a duration T+ and decreases from A+ for a duration T. The values for these four unknowns are found by enforcing periodicity of the solution and requiring it to cross-threshold twice, giving us four simultaneous equations:Embedded Image(4.11)Embedded Image(4.12)Embedded Image(4.13)Embedded Image(4.14)We solve these to findEmbedded Image(4.15)assuming 1>h (so that threshold can be reached). The amplitudes A± satisfyEmbedded Image(4.16)where T=T++T is the period of oscillation. We thus find that T satisfies the transcendental equationEmbedded Image(4.17)where R=1/h. The absolute amplitude A=A+A is given by Embedded Image. A plot of the period and amplitude as a function of τd is shown in figure 5. By linearizing about the periodic orbit shown in figure 4 and finding its Floquet exponents, one can show that this orbit is actually unstable (Coombes & Laing in press).

Figure 4

A periodic solution in a single population model with excitatory self-feedback. In this example, a=1, b=0, −θu=h=0.5 and τ1=τd=2.

Figure 5

The period and amplitude of an oscillatory solution in a single population with excitatory self-feedback as a function of the delay τd. The other parameters are the same as shown in figure 4.

5. Numerical bifurcation analysis

In the high-gain limit (when f is the Heaviside), explicit solutions of (2.1) can be constructed, as in §4. For a general firing rate function, solutions cannot normally be explicitly constructed, but the bifurcations of fixed points can be detected and followed in parameter space, as in §3. DDE-BIFTOOL (Engelborghs et al. 2001, 2002) is a software package for the numerical bifurcation analysis of systems of delay differential equations that can not only detect the bifurcations of fixed points, but can also follow branches of stable and unstable periodic orbits, and homoclinic and heteroclinic orbits. In this section, we demonstrate its capabilities by analysing (2.1) as θu and τ1=τ2τ are varied. Typical results are shown in figure 6, where the curves of saddle-node and Hopf bifurcations of fixed points are shown, along with the saddle-node bifurcations of periodic orbits. Here, as expected from §3, varying τ does not change the fixed points, but it does affect their stability. Figure 7 shows horizontal slices through figure 6 at τ=0.5, 0.2 and 0.09. For τ=0.5, there is a branch of stable periodic orbits joining the Hopf bifurcations on the upper and lower branches of fixed points. Between τ=0.5 and 0.2, a pair of saddle-node bifurcations of periodic orbits is created, resulting in the creation of a branch of unstable periodic orbits. For τ=0.09, an unstable periodic orbit is created from the Hopf bifurcations on the unstable middle branch of fixed points.

Figure 6

Bifurcation diagram. Solid line, saddle-node bifurcation of fixed points; dashed line, Hopf bifurcation; circles joined by a line, saddle-node bifurcation of periodic orbits. The parameter values are α=1, θv=0.5, τ1=τ2=τ, β=60, a=−1, b=−0.4, c=−1 and d=0.

Figure 7

Horizontal cuts through figure 6 at (a) τ=0.5, (b) τ=0.2 and (c) τ=0.09. Solid line, stable fixed points; dashed line, unstable fixed points; circles, stable periodic orbit; crosses, unstable periodic orbit (the maximum of u over one oscillation is plotted). The parameter values are the same as shown in figure 6. Note the different axis scales.

Brute force numerical simulation can also be used to explore small systems of delay differential equations. For example, Battaglia et al. (2007) studied a system very similar to ours, setting a=d<0 and b=c>0, but using a threshold linear firing rate function: f(z)=z if z>0, and zero otherwise. They varied both local and long-range interaction strengths (a and b in our notation) and found various types of chaotic and periodic behaviour. We have performed a similar calculation, with results shown in figure 8. For these parameter values, the system appears to have only one fixed point, and this undergoes a Hopf bifurcation on the curve shown. The most positive Lyapunov exponent can be found in the same way as for a system of ordinary differential equations (ODEs), by numerically integrating the variational equation in parallel with the underlying system. (Note that only one initial condition was used for each point in the parameter space, so multistability is not detected.)

Figure 8

Maximal Lyapunov exponent. The black line marks a Hopf bifurcation, to the right of which there is a stable steady state. A positive exponent indicates chaotic behaviour. The parameter values are α=1, θu=θv=0.2, τ1=τ2=0.1, β=60, a=d and b=c.

Figure 9a shows a typical chaotic solution corresponding to the point (a, b)=(−6, 2.5) in figure 8. Figure 9b shows a quasi-periodic orbit that was obtained using the parameter values in figure 9a, but simply decreasing β (the steepness of the firing rate function).

Figure 9

(a) A chaotic solution and (b) a quasi-periodic solution. The parameters are α=1, a=d=−6, b=c=2.5, θu=θv=0.2 and τ1=τ2=0.1, with (a) β=60 and (b) β=40.

6. Discussion

The periodic and chaotic behaviours of the type seen above are of great interest in neural systems, as are ‘bursting’ oscillations (Coombes & Bressloff 2005). Although the origins of bursting in low-dimensional ODEs are quite well understood, there has been very little work on bursting in delay differential equations. Here, we briefly summarize the results of several groups. Destexhe & Gaspard (1993) studied a system of two coupled DDEs, meant to model interacting populations of excitatory and inhibitory neurons. By varying one parameter they found bursts containing different numbers of APs. The bursting could be understood as resulting from a homoclinic tangency to an unstable limit cycle, and did not require the usual ‘slow–fast’ analysis (Coombes & Bressloff 2005). When the delays in their system were set to zero, the bursting could not exist, since the system was then two-dimensional. However, the general presence of a delay is not necessary to observe this bifurcation, as it can appear in three-dimensional ODEs (Hirschberg & Laing 1995).

Laing & Longtin (2003) studied the effects of paired delayed excitatory and inhibitory feedback on a single integrate-and-fire neuron, with and without noise. By assuming that the feedback was slow relative to the membrane time constant, they derived a rate model for the dynamics. With either inhibitory or paired excitatory and inhibitory feedback, these authors found periodic and chaotic oscillations in the firing rate of the neuron, i.e. bursting. They verified many of their results by simulating an actual integrate-and-fire neuron with appropriate delayed feedback.

Throughout this paper we have focused on discrete delays in neural population models without spatial extent. However, there is a large body of literature devoted to continuum models of neural tissues, particularly with regard to understanding the mechanisms of pattern and wave formation (see Coombes (2005) for a review). Many of the techniques we have touched upon here may be adapted for the treatment of such neural field equations (which are typically written as non-local evolution equations of integral type). Indeed, work in this direction has already been pursued by Roxin et al. (2005) in the context of macroscopic pattern formation in the cortex, and by Golomb & Ermentrout (1999) and Bressloff (2000) for the analysis of travelling waves in synaptic networks of integrate-and-fire neurons. More recent work on space-dependent delays (induced by the finite conduction speeds of APs along axons) can be found in Atay & Hutt (2006), Laing & Coombes (2006) and Coombes et al. (2007).

In summary, delays are ubiquitous in neural systems and should therefore be included in any realistic neural model. Here, we have briefly outlined the types of analysis available for small systems of neuronally inspired delay differential equations. There remains much to be discovered about the role of delays in more realistic neural models.


  • One contribution of 10 to a Theme Issue ‘Delay effects in brain dynamics’.

  • This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


View Abstract