## Abstract

We use weakly coupled oscillator theory to study the effects of delays on coupled systems of neuronal oscillators. We explore, first, simple pairs with constant delays and then examine the role of distributed delays as would occur in systems with dendritic branches or in networks where there is a distance-dependent conductance delay. In the latter, we use mean field theory to show the emergence of travelling waves and the loss of synchronization. Next, we consider phase models with stronger coupling and delays in the state variables. We show that they have a richer dynamics but one that is still similar to the weakly coupled case.

## 1. Introduction

Rhythms are ubiquitous in the nervous system and range from at least 200 cycles s^{−1} to less than one cycle per day (Yuste *et al*. 2005; Buszaki 2006). These rhythms are generated by diverse mechanisms such as the fast voltage-gated ion channels in the nerve membrane or the very slow molecular interactions that generate, say, the 24 hour circadian rhythm. Most of these rhythms (at least in vertebrates) are generated by the coherent oscillations of many coupled neurons. Sometimes they are synchronous and in other cases they are organized into patterns such as waves (for example in the swimming pattern of the lamprey; Cohen *et al*. 1992) or even rotating waves (in certain brain slice preparations; Huang *et al*. 2004.) Thus, one of the outstanding theoretical questions in our understanding of neural rhythms is how the coupling interacts with the intrinsic dynamics to generate these patterns of oscillations. The methods for studying interactions between coupling and intrinsic dynamics and what their consequences are for rhythmic behaviour invoke a variety of tools from nonlinear dynamics and applied mathematics. Two of the most commonly used methods are singular perturbation and the method of averaging (see Borisyuk *et al*. 2005, chs. 2 and 3). The latter is the focus of this paper.

Synchronous oscillations at frequencies ranging between 5 and 80 Hz are ubiquitous in the brain (Buszaki 2006). The production and synchronization of these rhythms are an area of active research with models ranging from the very complex (Traub *et al*. 2005) to more abstract (Ermentrout & Kleinfeld 2001; Izhikevich 2007). In order for any rhythms to synchronize, they must either be somehow connected or driven by common input. In order to synchronize a pair or more of connected oscillators, information about one oscillator must be communicated to the other. In the nervous system, the dominant means of communication is via synaptic transmission in which the firing of one neuron results in a perturbation of the voltage of the other neuron. The voltage perturbation interacts with the dynamics of the oscillator to shift it and, if the shift is in the appropriate direction, the neurons will tend towards synchrony. However, most neuronal communication is not instantaneous; rather, there are delays in the interactions. There are many sources of delays in the nervous system, ranging from conductance delays as action potentials travel down nerve axons to delays as neurotransmitter crosses the synaptic cleft. Other less obvious delays are found in the transmission of synaptic inputs through the long dendritic structures found in many pyramidal and other CNS neurons. Figure 1 illustrates some of these sources in a ‘typical’ pyramidal cell (named for its shape) from the cerebral cortex. The information flow is from A to B to C to A. In figure 1A, the signal from a different neuron produces a change in the membrane potential, *V*, at the site of the connection (synapse) depicted by A1. This potential is felt at the soma of the cell as a filtered (temporally shifted and attenuated) change in the potential (A1–A3). Thus, depending on where the input comes into the dendrite, there can be quite a substantial change in the shape of the interaction. Once the neuron generates an action potential or spike, it must send this down a long axonal fibre (figure 1B1–B3). Since this spike propagates at a finite conduction speed (of the order of 1 m s^{−1}), there can be a lengthy delay for long axons. For example, transmission of a spike from one hemisphere of the brain to the other can take as long as 5–10 ms (Whittington *et al*. 2000). For a 40 Hz rhythm, this is a substantial part of the cycle. Once the action potential reaches its target (figure 1C), it causes a release of neurotransmitter that must diffuse across the synaptic cleft and bind to the dendrite of the other neuron. The time from release to binding is yet another source of delay. We can view these three sources of delay as (i) distributed delay owing to the properties of dendrites, (ii) distance-dependent delays owing to finite propagation speed, and (iii) fixed delays owing to release and diffusion.

The general analysis of systems of coupled nonlinear delay differential equations is a formidable task. One reason for this is that they represent infinite-dimensional dynamical systems. However, there is at least one context in which this complexity can be significantly reduced. When oscillators are coupled together with weak interactions, then it is possible to use the theory of averaging to reduce the full coupled model to a set of equations on a torus (Ermentrout & Kopell 1984; Izhikevich 1998). In this scenario, delayed interactions between the oscillators are manifested as simple phase shifts in the interactions on the torus (see below). When delays are of the order of 1/coupling-strength, then it is no longer possible to reduce them to phase-shifts of the interactions and, instead, the delays appear inside the state variables (Izhikevich 1998). The goal of the rest of this paper will be to review some of the work that has been done in these two scenarios (‘regular’ and ‘large’ delays) both of which concern differential equations on a torus (also called phase models). We first set up the problem of neural communication in a fairly general setting and use the theory of averaging to reduce networks of coupled oscillators to networks of scalar phase equations. Then, we explore the consequences of the three types of delays described in the previous paragraph, first for regular delays and then for large delays. We will make the meaning of regular and large precise in §2.

A number of people have studied the role of delays in both the desynchronization and synchronization of oscillators (Schuster & Wagner 1989; Niebur *et al*. 1991; Nakamura *et al*. 1994; Bressloff & Coombes 1997*a*,*b*; Crook *et al*. 1997; Kim *et al*. 1997; Yeung & Strogatz 1999; Choi *et al*. 2000; Zanette 2000; Jeong *et al*. 2002; Earl & Strogatz 2003; Ko *et al*. 2004; Ko & Ermentrout 2007). Many of these papers use phase models with explicit delays (large case) which arise when the delays are long in comparison to the coupling strength. Schuster & Wagner (1989) show multiple stable synchronous and anti-phase states, while Niebur *et al*. (1991) show disorganized states in two-dimensional nearest neighbour coupled arrays. Asynchronous states in all:all coupled systems with delays are studied in Kim *et al*. (1997), Yeung & Strogatz (1999) and Choi *et al*. (2000). Bressloff & Coombes (1997*a*,*b*), Crook *et al*. (1997), Zanette (2000), Ko *et al*. (2004) and Ko & Ermentrout (2007) show that waves can be formed in the presence of distance-dependent delays. Earl & Strogatz (2003) provides a general condition for the stability of synchrony in the presence of explicit delays. Here, we review some of this work, provide some new examples and describe the general theory of weak coupling and synchrony in the presence of various types of delay.

## 2. Weak coupling

The modelling of individual neurons ranges from treating each cell as a simple point to including all the complex structure of the neuron. In this paper, we will primarily focus on point neurons to derive phase models and then extend these to models that have a very simple spatial structure such as a single linear dendrite. We will start with a derivation for a pair of coupled point neurons. Recall that models of neurons involve current balance equations corresponding to the many different ions involved in creating the action potential, the input currents owing to other neurons, and capacitative currents owing to the nature of the membrane (see Dayan & Abbott 2001, ch. 5 for an introduction). Consider a pair of neural oscillators that are coupled via a chemical synapse(2.1)Here, *V*_{i} is the voltage of the neuron and *s*_{j}(*t*) is some function or functional of the neuron labelled *j* (called the presynaptic cell, while the neuron *i* is called the postsynaptic cell). *C* is the membrane capacitance that we assume to be 1 without loss of generality. The function *I*_{ionic} represents all the currents that cause the neuron to spike rhythmically and the terms delineated by dots represent all the gating variables, ionic species, etc. *g*_{ij} is the coupling strength and *I*_{i} represents possible heterogeneities in the oscillators (e.g. applied currents). (Here, as in many biophysical models, voltage is measured in millivolts, currents in μA cm^{−2}, conductances in mS cm^{−2}, time in milliseconds and capacitance in μF cm^{−2}.) Setting *I*_{i}, *g*_{ij} to zero, we assume that the voltage and other variables (including the functional, *s*) have an asymptotically stable limit cycle with period, *T*. With no coupling, the pair of oscillators can be described by a pair of variables (*θ*_{i}, *θ*_{j}) each lying on a circle of circumference *T*. We set *θ*=0 as the time of the spike of the neuron. In the absence of coupling, the dynamics lies on an attracting torus. Thus, for sufficiently small coupling, the torus will persist and our interest is in the dynamics on this torus. It can be readily shown (Ermentrout & Kopell 1984) that the dynamics satisfies(2.2)where(2.3)Here, *Δ*(*t*) is the phase resetting curve for the neuron subjected to a brief current pulse. (More correctly, *Δ*(*t*) is the voltage component of the periodic solution to the adjoint of the linearized equation for the neuron evaluated along the limit cycle (Ermentrout & Kopell 1984; Izhikevich 2007). Suffice it to say that it characterizes how the oscillator is phase shifted due to small perturbations at time *t* after the last spike.) The constant *K* is just the mean value of *Δ*(*t*). Equation (2.3) is the key to our understanding of how delays and other linear functionals alter the coupling between neurons. Let us for the moment drop the indices in the function *H*_{ij} since if the two neurons are otherwise identical, then in the absence of coupling, *V*_{i}(*t*)=*V*_{j}(*t*), *s*_{i}(*t*)=*s*_{j}(*t*) up to a phase shift. Let us suppose that *s*(*t*) is some filtered or delayed version of a function *S*(*t*) (that we can call the instantaneous interaction—the one that occurs if there are no delays of any kind) and define (we use *R*(*t*) to save on notation)We next define the effective interactionSuppose we can writeFor example, if *k*(*t*)=*δ*(*t*−*τ*), then *s*(*t*)=*S*(*t*−*τ*), that is, *s* is a delayed version of *S*. Using the definition of *H*_{s}, we see thatWe interchange the order of integration and observe that(2.4)That is, any linear functional of the instantaneous synaptic input becomes the same linear functional of the associated interaction function. For example, in the case of a delay, *H*_{s}(*ϕ*)=*H*_{S}(*ϕ*−*τ*). All that the delay does is shift the interaction function.

In a network of oscillators with delays *τ*_{ij} between oscillator *j* and oscillator *i*, and applied currents *I*_{i} to account for heterogeneities, we thus have the following system of general phase equations:Here, we have assumed that other than the strength and delay, all interactions are the same. To emphasize the smallness of the interactions and heterogeneities, we write and where the are *O*(1) and 0<*ϵ*≪1. Finally, we let where is a slow time scale. With this notation, the reduced equations become(2.5)Henceforth, we will drop the hats from all the variables and parameters. This derivation is good to order *ϵ*.

For fixed constant delays, Izhikevich (1998) also derives this equation but shows that there is an additional term that could be of the same order of magnitude when the delays are large. He shows that equation (2.5) can be written (without the hats) as(2.6)Note that there is an explicit delay of *θ*_{j}(*t*) by an amount *ϵτ*_{ij}. If *τ*_{ij} is *O*(1), then this term contributes very little to the interactions and can be ignored. However, if the delays are large compared to the coupling, then *ϵτ*_{ij} could be *O*(1), and thus contribute to the dynamics. Our goal in the rest of this paper is to review the work that has been done in equations (2.5) and (2.6), which we define to be the regular and large delay cases, respectively. In order to formally be in the large delay case, the delays must be of the order of 1/*ϵ*; for weak coupling, then, the delays must be quite large. In spite of this restriction, there are many papers on phase models with explicit delays that are mathematically interesting in their own right. However, if you want to connect equations of the form (2.6) to a biophysical model of the neuron, you must realize that you are formally assuming very large delays.

## 3. Regular delays

In this section, we explore the effects of *O*(1) delays on the weak coupling theory. We first consider a pair of oscillators with constant delays and then the same pair when the delay is the result of passive filtering through dendrites. Then, we look at a chain of oscillators with random coupling but distance-dependent delays and show the existence of wave-like behaviour.

### (a) Pairs with constant delays

Let us consider the simplest possible network of two coupled oscillators that are identical and symmetrically coupled but have a fixed delay, *τ*, in their interaction:Let *ϕ*=*θ*_{2}−*θ*_{1} denote the phase difference between the two oscillators. Then(3.1)Zeros of *G*(*ϕ*, *τ*) represent phase-locked patterns for the pair of cells. Note that *G*(*ϕ*, *τ*) is proportional to twice the odd part of the *τ*-shifted interaction function *H* and so it is always an odd periodic function. This means that there are at least two roots, the synchronous solution, *ϕ*=0, and the anti-phase solution, *ϕ*=*T*/2. Synchrony is asymptotically stable if *G*_{ϕ}(0,*τ*)<0 which implies that *H*′(−*τ*)>0. This tells us an interesting result right away. Suppose that in the absence of delays, synchrony is stable; then, *H*′(0)>0. However, since *H* is periodic and (generally) continuously differentiable, there must be a point at which *H*′(−*τ*) becomes negative. This implies that for sufficient delay relative to the period of the oscillator, synchrony will be destabilized. A ‘rule of thumb’ which is often heard is that if the delay exceeds a quarter of the period, synchrony will be unstable. This, is not really true, but for *H*(*ϕ*)=sin(*ϕ*) (a commonly used interaction function), it is clearly the case that for delays between *π*/2 and 3*π*/2, synchrony will be unstable. Of course, the pure sine interaction function is an idealization, so we will focus on a ‘real’ interaction function. The Morris–Lecar (ML) model is a simple two-variable model for a spiking neuron which we can synaptically couple with explicit delays. For weak coupling, the interaction function is readily computed and we can test the role of delays in the synchrony of a pair. Equations for the model are in appendix A. Figure 2*a* shows *G*(*ϕ*, *τ*) for three values of the delay *τ* (measured in ms) in the ML model with a period of 91.25 ms. For *τ*=0, it is clear that the slope at the origin is negative and that at a half-cycle is positive so that, with no delay, synchrony is stable. When *τ*=20, synchrony is unstable (slope is positive) while the anti-phase solution is stable. Interestingly, when *τ*=10, there is bistability between synchrony and anti-phase. Figure 2*b* shows a (*V*_{1}(*t*), *V*_{2}(*t*)) phase plane (units are millivolts) for the full ML model coupled weakly (*g*=0.001 mS cm^{−2}) with a synaptic delay of 15 ms. Initial data are started near synchrony for the black diagonal curve and near anti-phase for the grey curve. Solutions converge to synchrony or anti-phase, respectively. Thus, there is bistability. By adding a bit of noise to the full equations, we can record the time that neural oscillator 2 fires relative to the last spike from oscillator 1. A histogram of these times is shown in figure 2*c* for three different delays. While figure 2*a* shows bistability for *τ*=10, the histogram shows a single broad peak at the anti-phase time. This is because the anti-phase state has a sharper slope than the synchronous state leading to a much greater probability of being near to the anti-phase state. For *τ*=8, the function *G* is quite small in magnitude so that the distribution of times is closer to uniform.

### (b) Dendritic delays

Many neurons have long dendritic processes through which the synaptic input is filtered (Bressloff & Coombes 1997*a*). Crook *et al*. (1998) were among the first to look at this; later Lewis & Rinzel (2003) applied the ideas to gap junctions between neurons. We sketch a slightly simpler scenario here in which we suppose that there is an infinite dendrite that has a small soma coupled at *x*=0, where the dendrite influences the somatic dynamics but such that the soma has no influence on the dendrite (this is not completely realistic, but if the size of the soma is small compared to the dendrite, the effects of the soma on the dendrite will be small). At a point *X*, there is a synapse from another cell onto the dendrite and the potential at the origin provides the coupling between the two neural oscillators. We assume that the synapse is in the form of a current rather than an increase in conductance (as this makes the analysis simpler; figure 3*a*). With no loss in generality, we translate the resting potential of the dendrite to 0 mV. Then the dendrite satisfies the differential equationwhere *I*_{pre}(*t*) is the current from the presynaptic cell. At the soma (here, the post-synaptic cell), the equations areHere, *g*_{DS} is the coupling conductance between the dendrite and the soma; this is small in order to exploit the weak-coupling assumption. To understand the coupling, we need to solve for the dendritic potential, *V*(*x*, *t*). We use the Green's function for the infinite dendrite (Dayan & Abbott 2001)where *C*_{0} is a constant that depends on the diameter of the dendrite and the transmembrane resistance. From this expression, we obtain *V*(*x*, *t*) at the soma asThe weak coupling expression for the interaction function is thusThe phase-dependent (*ϕ*-dependent) part of this is exactly as we wrote in equation (2.4)(3.2)Thus, the dendrite acts as a generalized filter on the interactions. Figure 3*b* shows the Green's function as a function of time at three different locations, *X*/*λ*=1, 2 and 4 for a dendrite with a time constant, *τ*_{m}=1. To clarify the shapes, we have magnified the curve so that they all have the same peak size. Figure 3*c* shows the effect of this distributed delay on a sinusoidal interaction function. Recalling that, for example, synchrony is stable if the slope of the interaction function is positive, we see that the position of the synapse on dendrite can change whether or not synchrony is a stable state for identical pairs. Applying this to the Morris–Lecar model, we can plot the slope of the effective interaction function (3.2) at synchrony (*ϕ*=0) and at the anti-phase solution (*ϕ*=*T*/2) as the position of the synapse moves up the dendrite away from the soma. This is done in figure 3*d*. Note that for very close synapses, synchrony is stable (the slope, *H*′(0), is positive, but anti-phase is unstable (*H*′(*T*/2)<0)). However, for synapses farther out, the situation is reversed and for a small region (approx. *λ*/2 away), neither synchrony nor anti-phase is stable and instead there is a non-zero stable phase difference between the two cells.

In summary, for a pair of mutually coupled neurons, delays can arise either directly or from the passive interactions at a distance.

### (c) Distance-dependent delays

Large randomly coupled networks of oscillators which would pairwise synchronize have a tendency to generate global synchronous activity since the average coupling felt by each oscillator is roughly the same. That is, if the coupling is 0 or 1/*N* with probability *p*, then, for *N* large, the average coupling felt by each oscillator is *p* with variance *O*(1/*N*). However, suppose we impose some geometry in the nature of the interactions and assume the same random coupling but such that there is a distance-dependent delay. For example, we can arrange the oscillators on a ring and make the coupling between them depend on the circular distance on the ring. Then, as we saw above, with sufficient delays, synchrony is disrupted and new kinds of patterns can emerge. Before considering these, we briefly discuss a different scenario. Suppose we start with a ring of oscillators such that each one interacts with 2*m* nearest neighbours. If *m*/*N* is small enough, and interactions are synchronizing, then, in addition to synchronous solutions, there are also a variety of stable travelling waves corresponding to one or more complete cycles around the ring (i.e. *θ*_{i}(*t*)=*Ωt*+2*π*i*k*/*N*, where *k*=0,±1, …). Starting with this kind of network, ‘small-world’ networks are generated (Watts & Strogatz 1998) by picking a connection between a pair and randomly reassigning it. For example, with probability *q*≪1, the connection that oscillator 3 receives from oscillator 4 is replaced by one from oscillator *j* where *j* is randomly chosen. Barahona & Pecora (2002) and others have shown that manipulating the network in this fashion greatly fosters the appearance of synchrony and drastically reduces the basin of attraction for wave-like behaviour. However, if the interactions between oscillators have finite communication speed, then a link at a distance may incur a large delay and thus the new long-distance pairwise interaction may be desynchronizing. Instead of encouraging synchronization, it might actually foster waves. Ermentrout & Kopell (1994), for example, showed that in a network of nearest neighbour-connected oscillators, the presence of one distant desynchronizing link can disrupt the synchronous solution and lead to patterned solutions. To illustrate this idea, we consider a ring of twenty oscillators with the nearest neighbour couplingSynchrony is an asymptotically stable state for this system. In addition, there is also a stable travelling wave solution, *θ*_{i}=*ωt*+2*πi*/20. Now, suppose we replace the connection that oscillator 1 receives from oscillator 20 with a connection from oscillator 1+*l*. If there is no delay due to the distance of this connection, then synchrony will remain a stable solution. However, as soon as this connecting incurs a delay, synchrony will not even be a solution and some other pattern will arise. Figure 4*a* shows an example of the disruptive effects of the delay on synchrony. The relative phases of all 20 oscillators are shown when the coupling from oscillator 20 to oscillator 1 is replaced by an oscillator *l* unit away that has a delay of 2. Curiously, the more local connections are more disruptive in the sense that the peak of the relative phases is bigger when *l*=2 than when *l*=10. Figure 4*b* shows the relative phase of oscillator 10 as the delay increases. For 0 delay, all oscillators are synchronized. As the delay increases, there is an almost linear increase in the relative phase with local disruptions *l*=2 having a greater slope than longer-range disruptions, *l*=4, 10. In all three cases, the patterns disappear at a saddle-node bifurcation (marked with an asterisk) as a particular delay. The unstable branch of solutions meets with another stable branch at higher delay (marked by a circle) and for each case there is a regime of bistability which, for larger *l* extends over the full range of delays. Figure 4*c* shows the consequence of switching connections for a travelling wave solution (the diagonal line in the figure). Here, for three of the curves, we do not delay the interaction; the wave becomes more perturbed, for larger values of *l*. However, adding the delay of 2 to the *l*=10 case makes the solution closer to the unperturbed wave (yellow curve: *l*=10, *τ*=2). Thus, delayed long-range connections are conducive to both disrupting synchrony and promoting waves even at the level of changing a single connection. We will see this more systematically in a sparsely connected network with distance-dependent delays.

Let us consider a random network with distance-dependent delays(3.3)Here is the mean number of connections to each oscillator. We have written this in such a way that we do not have to worry about the circular distance and we treat indices modulo *N*, so that *i*±*N* is equivalent to *i*. *τ*′ characterizes the degree of the delay. Furthermore, we scale the distance-dependent delay, by 1/*N*. Figure 5*a* shows a simulation of equation (3.3) where *N*=1600, *H*(*ϕ*)=sin *ϕ* and the connection probability between oscillators is *q*=0.025; so for 1600 oscillators, each is connected to roughly 40 other oscillators. Initial data consist of random assignment of phases between −*π* and *π*. For low delays (*τ*′ small), solutions go to synchrony (figure 5*a*(i)), while at higher delays, there are near travelling-wave solutions (figure 5*a*(i–iii)). To understand this, Ko & Ermentrout (2007) show that if *A*_{ij} is 1 with probability, *q*, then as *N*→∞ the behaviour of this model approaches that of the mean-field equations(3.4)Crook *et al*. (1997) considered an equation similar to (3.4) on the infinite line in a slightly different context. It is easy to see that there is a family of travelling wave solutions of the form *θ*(*x*, *t*)=*Ω*_{m}*t*+2*πmx*, for *m*=0, ±1, … withNote that one can think of the dependence of *Ω*_{m} on *m* as the dispersion relation for the medium. If *H* is an odd periodic function and *τ*′=0, then *Ω*_{m}=*ω* and the frequency of all waves is the same. What is interesting is that the delay can have a profound effect on the stability of waves and synchrony. The eigenvalues of the linearization of equation (3.4) about the *m* travelling wave, *λ*_{n,m}, satisfywhere *n* is an integer and *m* is the wavenumber of the travelling wave. Figure 5*b* shows the maximum over *n* of the real parts of these eigenvalues. Regimes where this is negative are stable. These integrals can be evaluated by writing *H* in a Fourier series and then evaluating the resulting integrals. Another way to establish stability is the long-wave approximation. For *n* large the cosine part of the integral is *O*(1/*n*), so that a necessary condition for stability is thatFigure 5*c* shows that this approximation is a good match for the full stability criterion shown in figure 5*b*. The agreement is qualitatively good and brings home the main point that the regions of stability alternate periodically.

In this section, we have shown that delays in weakly coupled arrays of oscillators can have profound effects in both the patterns of behaviour and in their stability. In §4, we explore the consequences of large delays, that is when the term *ϵτ*_{d} is *O*(1) in equation (2.6).

## 4. Large delays

So far, delays in the phase model have been implemented as mere phase shifts in the interaction functions as dictated by the theory of averaging. However, as Izhikevich (1998) shows, delays that are large enough should be included as explicit delays of the phase variable. In this section, we consider some example networks with explicit delays and we show that many of the results of the previous part of the paper still hold. We start with a simple example to show the complexities of the problem. Schuster & Wagner (1989) consider the following model:In the normal delay regime, we noted above that the *τ* does not explicitly appear in the phases *θ*_{j}(*t*), but rather is implemented as a shift in the sine function. Nevertheless, the results outlined above suggest the kind of solutions one could expect and we can use this to reduce the existence of phase-locked solutions to an algebraic problem. We suppose that *θ*_{1}(*t*)=*Ωt* and *θ*_{2}(*t*)=*Ωt*+*ϕ* where *Ω* is the ensemble frequency and *ϕ* is the phase difference between the two oscillators. We substitute this ansatz into the equations and find(4.1)(4.2)which is two equations in two unknowns. Before examining this equation, let us recall the regular delay case equationsThere are still two equations and two unknowns, but critically, the unknown *Ω* does not occur inside the sine function, so we can subtract these equations to get a single scalar equation for *ϕ*, *ω*_{2}−*ω*_{1}=2*K* sin *ϕ* cos *τ* that can then be substituted back to get a value of *Ω*. In particular, if *ω*_{2}=*ω*_{1}, the only phase-locked solutions are *ϕ*=0, *π* corresponding to synchrony and anti-phase. However, with a large delay, we cannot separate the two equations and thus must solve two nonlinear equations simultaneously. Schuster & Wagner (1989) provide a fairly complete numerical analysis of the problem including stability. We consider a particular case in which the oscillators are identical, *ω*_{1}=*ω*_{2}=*ω*. Then subtracting equation (4.1) from (4.2), we get after simplificationAs with the regular delay case, there are always the solutions *ϕ*=0, *π*. From this, *Ω* is determined from the equation(4.3)which has many solutions if *K* is large. In addition to these solutions, we can also choose *Ω*=*Ω*_{n} where *τΩ*_{n}=(2*n*+1)*π*/2 for *n* an integer and *ϕ* such thatFor each *K*, there are only finitely many roots since |*Ω*_{n}| increases linearly with *n* and thus we will not be able to find a root *ϕ*. Unlike the regular delay case, the number of phase-locked solutions depends strongly on the coupling strength. With large *K*, there are many different roots and a large degree of multistability.

We can easily analyse the linear stability of these solutions. The linearization isRecall that there are two sorts of solutions: those in which *ϕ*=0, *π* and those in which *τΩ*=(2*n*+1)*π*/2. In the latter case, there is always one positive eigenvalue, *λ*, satisfyingwhere *a*=sin *ϕ*. In the former case (*ϕ*=0, *π*), solutions are stable as long as cos(*ϕ*−*Ωτ*)>0. Note that in the absence of delays, *ϕ*=0 is always stable. If the coupling *K* is small, then *Ω*≈*ω* and we can see that for the usual value of *ω*=1, the stability condition is the same as we derived in §3*a*. These results show that for large *K* there can be multistability between fast and slow synchronous and anti-phase solutions.

More generally, Earl & Strogatz (2003) consider the following system of delay equations:where *a*_{ij} is 1 for *k* elements and 0 otherwise. There is a synchronous solution, *θ*_{i}(*t*)=*Ωt*, provided(4.4)The left-hand side is a straight line and the right is a periodic function whose amplitude is dependent on the coupling strength. Thus, if *K* is very large, there will be many intersections and thus many possible frequencies, *Ω*. Earl & Strogatz (2003) show that such a synchronous solution is stable if and only if *KH*′(−*Ωτ*)>0. This condition is the same as found for the simple pairwise sinusoidal model discussed above. Thus, in general when coupling is strong there may be multiple stable synchronous solutions.

For general coupling of general phase model oscillators with delays, one can look for phase-locked solutions of the form *θ*_{i}=*Ωt*+*ϕ*_{i} where *ϕ*_{1}=0 and *Ω*, *ϕ*_{2}, …, *ϕ*_{N} are constants. For example, Ko *et al*. (2004) considered the large-delay analogue of equation (3.3)where *r*_{ij} is the distance between two coupled cells; *v* is conduction velocity; *A*_{ij} is 0 or 1 according to whether the oscillators are connected; and *n*_{i} is the number of inputs into oscillator *i*. The geometry is a circle of length *L* with each oscillator placed at *Li*/*N*. As in figure 5 above, they find the existence of stable travelling waves, *θ*_{i}≈*Ωt*+2*πmi*/*N*, where *Ω* approximately satisfiesNote that *Ω* appears on both sides of the equation, so that unlike the situation discussed above, computation of *Ω* is difficult.

Jeong *et al*. (2002) consider other geometries with delay including two dimensions where they find patterned states that include plane waves, target waves and spiral waves. Specifically, they consider an array of *N*×*N* oscillators coupled in a periodic domain with interactions a given radius(4.5)Here *ω* is the natural frequency of oscillators; *K* is the coupling strength; ; and , where *r*_{kl,ij} denotes the Euclidean distance between (*i*, *j*) and (*k*, *l*): . *v* is the velocity of transmission between the oscillators so that the delay is distance-dependent. Figure 6*b* shows a snapshot in time for *K*=1, *r*_{0}=4, *v*=1. This means that oscillators which are, say, 4 away will be delayed by four time units. In this particular case, spiral waves arise from random initial data but the spirals do not have a fixed axis of rotation. Rather, they drift, expand and contract. This pattern of behaviour as well as other simpler patterns can be easily understood by looking at a simpler model with the nearest neighbours and a fixed delay(4.6)where (*i*′, *j*′) are the four nearest neighbours of (*i*, *j*). Here, as in §3, the delay appears only as a phase-shift. Figure 6*a*,*c* show solutions to (4.6) for zero and periodic boundary conditions, respectively, as *τ* increases from 0. Paullet & Ermentrout (1994) proved that for the zero boundary condition case and pure sinusoidal coupling, there is a stable rotating wave solution. This is shown in panel *a*(i). As *τ* increases (*a*(ii–iv)), the rotating wave develops a twist and looks as if a spiral. The squaring of the arms is a consequence of the zero boundary conditions. In each panel of figure 6*a*, the delay is *τ*=*π*(*j*−1)/20 which is sufficiently small such that the solutions are phase-locked. That is, the relative phases shown in the figure are steady states. Panels *c*(i–iv) show a similar simulation with periodic boundary conditions. Note that four singularities (spirals) develop owing to the periodic conditions. In *c*(v), the delay is *π*/4 that results in a complex meandering and expanding set of spirals just as in figure 6*b*. Ermentrout (1995) provides a heuristic mechanism for the tightness of the spiral and the eventual loss of locking as a function of *τ*, the shift in the interaction function.

## 5. Discussion

In this paper, we have reviewed the consequences of various types of delays on the behaviour of networks of coupled phase models. Phase models represent a particularly simple class of models for interacting nonlinear limit-cycle oscillators that exhibit richness in behaviour while at the same time admit analytic approaches. Their most attractive feature is that they can be directly linked to more complex biophysical models. Modest delays between oscillators do not affect the structure of phase models but can have dramatic effects on their behaviour as they shift the interactions. This shift can result in destabilization of states that would otherwise be stable and when incorporated into spatially distributed models can lead to the emergence of new stable patterns of activity such as waves.

One issue that remains to be resolved is whether or not the phase models with explicit delays have behaviour other than phase-locking. For example, consider the simplest possible modelWe have already studied phase-locked solutions of this system of the form *x*_{i}(*t*)=*Ωt*+*ϕ*_{i} and found that there can be many such solutions depending on *τ* and *K*. However, none of these are qualitatively different from solutions to the *τ*-shifted phase models. As far as we know, there are no other attractors to this system. If sin(*θ*) is replaced by a more complex function, *H*(*θ*), our extensive simulations seem to indicate that there are still only phase-locked solutions.

We have neither addressed delays in the individual dynamics of the oscillators, nor have we discussed delay-induced oscillations (Glass & Mackey 1988; Doiron *et al*. 2003). Indeed, as far as we know, weak coupling theory has not yet been developed or applied to coupled delay equations (as opposed to delayed coupling, which we review in this paper). Thus, there are a number of open problems that remain in the area of coupled oscillators with delays.

## Footnotes

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

- © 2009 The Royal Society