## Abstract

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 form(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 sigmoid(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 equations(2.3)where . The Jacobian matrix is therefore(2.4)Thus, the conditions for a Hopf bifurcation (HB) are(2.5)Eliminating *v*^{*} as(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.

## 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 . Demanding that the amplitudes be non-trivial gives a condition on *λ* that may be written in the form (*λ*)=0, where(3.1)and the equilibrium (*u*^{*}, *v*^{*}) is given by the simultaneous solution of (2.3). For *λ*∈, we see that *λ*=0 when(3.2)where *κ*_{1}=*aβu*^{*}(1−*u*^{*}), *κ*_{2}=*dβu*^{*}(1−*v*^{*}) and *κ*_{3}=*bcβ*^{2}*u*^{*}*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 *ω*∈. The bifurcation condition in this case is defined by the simultaneous solution of the equations Re (i*ω*)=0 and Im (i*ω*)=0, namely(3.3)(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 . 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).

## 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).

### (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 3*a*, we parametrize such a solution in terms of two fundamental times *T*_{1,2} and the maxima and minima *A*_{±} of the orbit. Here, *T*_{1} denotes the time spent on the decreasing part of the trajectory, and *T*_{2} that spent on the rising phase. Exploiting the piecewise linear nature of the dynamics, we then have that(4.1)(4.2)(4.3)(4.4)Solving these we obtain the period of oscillation *T*=*T*_{1}+*T*_{2}, where(4.5)and . 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 3*b*, we note that by symmetry, the rising and falling phases have the same duration, say *T*_{1}. For the parameters considered, we find that *τ*_{1}<*T*_{1}<*τ*_{2}, and we obtain the relationships(4.6)(4.7)(4.8)Solving the above we find that *T*_{1} satisfies the transcendental equation(4.9)The period *T* is 2*T*_{1} and the absolute amplitude of oscillation, *A*=*A*_{+}−*A*_{−}, is given by(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:(4.11)(4.12)(4.13)(4.14)We solve these to find(4.15)assuming 1>*h* (so that threshold can be reached). The amplitudes *A*_{±} satisfy(4.16)where *T*=*T*_{+}+*T*_{−} is the period of oscillation. We thus find that *T* satisfies the transcendental equation(4.17)where *R*=1/*h*. The absolute amplitude *A*=*A*_{+}−*A*_{−} is given by . 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).

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

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

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

## Footnotes

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.

- Copyright © 2009 The Royal Society