Delays in activity-based neural networks

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.


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.

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 Here, u and v represent the synaptic activity of the two populations, with a relative time scale for the response set by a K1 . The architecture of the network is fixed by the weights a, b, c, d, while q 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 0 Zbf(1Kf ), with bO0. The fixed delays t 1 and t 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, bZ0, (ii) Chen & Wu (1999), when aZ1, aZdZ0 and bZcO0, and (iii) Battaglia et al. (2007), when aZ1, aZd!0, bZcO0 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 t 1 Zt 2 Z0, it is straightforward to find the values of q u and q 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 Thus, the conditions for a Hopf bifurcation (HB) are Tr L ZKð1 C aÞ C abu Ã ð1Ku Ã Þ C adbv Ã ð1Kv Ã Þ Z 0 and det LO 0: ð2:5Þ we can then plot the fixed point equation (parametrically) in the (q u , q v ) plane, as in figure 1. A similar procedure can be used to determine the locus of saddlenode (SN) bifurcations defined by det LZ0, as well as the Bogdanov-Takens bifurcation defined by det LZ0 and Tr LZ0 (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.

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 ðu; vÞZ ð u; vÞe lt . Demanding that the amplitudes ð u; vÞ be non-trivial gives a condition on l that may be written in the form E(l)Z0, where and the equilibrium (u Ã , v Ã ) is given by the simultaneous solution of (2.3). For l2R, we see that lZ0 when Thus, a real instability of a fixed point is defined by (3.2) and is independent of (t 1 , t 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 lZiu for us0, where u2R. The bifurcation condition in this case is defined by the simultaneous solution of the equations Re E(iu)Z0 and Im E(iu)Z0, namely 0 Z ð1K k 1 cosðut 1 ÞÞð1K k 2 cosðut 1 ÞÞKðu C k 1 sinðut 1 ÞÞðu=a C k 2 sinðut 1 ÞÞ K k 3 cosð2ut 2 Þ; ð3:3Þ For the parameters that ensure us0, we shall say that the simultaneous solution of equations (3.3) and (3.4) defines a Hopf bifurcation at ðt 1 ; t 2 ÞZ ðt c 1 ; t c 2 Þ. More correctly, we should also ensure that as the delays pass through this critical point that the rate of change of Re l 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 (t 1 , t 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).

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 aZ1, q u Zq v , aZdZ0 and bZc 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 (b/N), so that f(z)ZH(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 q u Zq v and aZ1. Regarding a synchronous T-periodic solution, u(t)Zv(t) with u(tCT )Zu(t), like that shown in figure 3a, we parametrize such a solution in terms of two fundamental times T 1,2 and the maxima and minima A G 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Þ Solving these we obtain the period of oscillation TZT 1 CT 2 , where ð4:5Þ and sZKðae t 1 C be t 2 Þ. The amplitude of the oscillation is AZA C KA K Z (aCbCs)/s. Similarly, to analyse an anti-synchronous solution, u(t)Zv(tCT/2) with u(t)Zu(tCT ), as in figure 3b, we note that by symmetry, the rising and falling phases have the same duration, say T 1 . For the parameters considered, we find that t 1 !T 1 !t 2 , and we obtain the relationships ð4:6Þ Solving the above we find that T 1 satisfies the transcendental equation The period T is 2T 1 and the absolute amplitude of oscillation, AZA C KA K , is given by 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 aZ1, bZ0, q u ZKh, hO0 and t 1 Zt 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 G and T G , which denote the largest (A C ) and smallest (A K ) values of the trajectory and the times spent above (T C ) and below (T K ) the threshold h. The trajectory increases from A K for a duration T C and decreases from A C for a duration T K . 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: We solve these to find assuming 1Oh (so that threshold can be reached). The amplitudes A G satisfy  Figure 4. A periodic solution in a single population model with excitatory self-feedback. In this example, aZ1, bZ0, Kq u ZhZ0.5 and t 1 Zt d Z2.
where TZT C CT K is the period of oscillation. We thus find that T satisfies the transcendental equation where RZ1/h. The absolute amplitude AZA C KA K is given by AZ ½e ðTKt d Þ K 1. A plot of the period and amplitude as a function of t 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).

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(Engelborghs et al. , 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 q u and t 1 Zt 2 ht 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 t does not change the fixed points, but it does affect their stability. Figure 7 shows horizontal slices through figure 6 at tZ0.5, 0.2 and 0.09. For tZ0.5, there is a branch of stable periodic orbits joining the Hopf bifurcations on the upper and lower branches of fixed points. Between tZ0.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 tZ0.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 aZd!0 and bZcO0, but using a threshold linear firing rate function: f(z)Zz if zO0, 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 9a shows a typical chaotic solution corresponding to the point (a, b)Z (K6, 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 b (the steepness of the firing rate function).

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