## T Fig. 1. Illustration of the stability and bifurcation results for the example of (12). The equilibrium solution is stable in the region contiguous with the t axis. The number of eigenvalues with positive real part is shown in each subregion of the plane. Thick/thin curves correspond to Hopf bifurcations giving rise to synchronous/ anti-phase oscillation

Fig. 1. Illustration of the stability and bifurcation results for the example of (12). The equilibrium solution is stable in the region contiguous with the t axis. The number of eigenvalues with positive real part is shown in each subregion of the plane. Thick/thin curves correspond to Hopf bifurcations giving rise to synchronous/ anti-phase oscillation

Recall that the equilibrium points of (13) with t > 0 are the same as those with t = 0. Thus for the neural model (13) with t > 0, the bifurcations involving only equilibrium points (saddle-node, pitchfork, transcritical) will be the same as those for (13) with t = 0.

The two main bifurcations leading to the creation of periodic orbits in neural systems are the Hopf bifurcation and the infinite period bifurcation. These bifurcations are associated with Type II and Type I oscillators, respectively (Breakspear and Jirsa, 2006).

Consider first the Hopf bifurcation. This involves the creation of a periodic orbit as an equilibrium point changes stability. There are simple criteria to check to determine if a Hopf bifurcation occurs in a delay differential equation at a particular parameter value, say t = tc.

### Hopf Bifurcation Test

Assume that system (13) has an equilibrium point x. If the following are satisfied, then system (13) undergoes a Hopf bifurcation at x as t passes through tc.

1. The characteristic (19) of the linearization of (13) about x has a pair of pure imaginary eigenvalues, when t = tc, that is,

2. As t passes through tc the rate of change of the real part of this eigen value^) is nonzero, that is, dRdê(X)

3. The characteristic (19) of the linearization of (13) about xv has no other eigenvalues with zero real part.

Other than in some exceptional cases, this is enough to guarantee that a periodic orbit is created as t passes through Tc.

Whether the periodic orbit is stable or unstable depends on the nonlinear terms in the equation. There are two main approaches for determining this analytically, both of which require intensive computations and are best done either numerically or with a symbolic algebra package such as Maple. The centre manifold construction reduces the system of delay differential equations to a system of two ordinary differential equations from which the stability of the periodic orbit (for t close to Tc) may be deduced. See (Belair et al., 1996; Wischert et al., 1994; Wu et al., 1999) for examples of how this is done. Perturbation methods, such as averaging and the method of multiple scales, find an approximate expression for the periodic solution and for the corresponding Floquet exponents. See (Campbell et al., 2006; Gopalsamy and Leung, 1996; Wirkus and Rand, 2002) for examples of how this is done.

Example. Applying this test to our coupled Fitzhugh-Nagumo model shows that the system has a Hopf bifurcation along each of the curves where the characteristic equation has a pair of pure imaginary eigenvalues, i.e., along the curves defined by (21)-(22) and shown in Fig. 1. By analyzing the solutions of the linearization (16) that correspond to the roots, one can show that some of the Hopf bifurcations give rise to synchronous or in-phase oscillations (i.e. vi(t) = v2(t) and wi(t) = w2(t) for all t) and some to anti-phase solutions (i.e. the spikes in vi and v2 are half a period apart and similarly for wi and w2 ).

One important thing to note about Hopf bifurcation in systems of delay differential equations is that there are always multiple branches of Hopf bifurcation. This can be seen in our example. The t value where a Hopf bifurcation occurs corresponds to a t value satisfying (22). Clearly if a given value of t satisfies this equation, then so does t + kn, k = ±1, ±2,

Now consider the the infinite period bifurcation. This bifurcation occurs when a saddle-node bifurcation occurs on an invariant circle. As indicated above, the conditions for the saddle-node bifurcation to occur in a delay differential equation are the same as for the corresponding system with zero delay. Whether or not this bifurcation occurs on a limit cycle is not easily determined analytically (even without delays), thus these bifurcations are often investigated using numerical tools (see Sect. 2.5).

### 2.3 Lyapunov Theory

The basic idea of Lyapunov theory is to use an auxiliary function to determine the dynamics of a nonlinear system. A very simple example is the total energy in a mechanical system with damping, such as the pendulum model:

The total energy of this system is

A simple calculation, keeping in mind that 0 and 0 depend on t, shows that ddE < 0. This means that as t increases, E must tend to a minimum value. This in turn determines what the solutions of the nonlinear model can do. In particular, one can show that this implies that all solutions must tend to one of the equilibrium points (0, 0) = (2kn, 0), k G ZZ as t ^ to, i.e. the pendulum swings with smaller and smaller amplitude until it is hanging straight down. Lyapunov theory generalizes this idea to an arbitrary auxiliary function, V (x), which has similar properties to the energy function in the above example, viz.,

1. V(x) > 0, x = 0; V(0) = 0 (V positive definite)

These properties can be used to show that the equilibrium point x = 0 is asymptotically stable. By modifying the properties above, one can also use Lyapunov functions to show that an equilibrium point is unstable, that all solutions are bounded or that all solutions synchronize as t ^to.

There are two ways of extending the Lyapunov theory for ordinary differential equations to delay differential equations such as (13). Lyapunov functionals are auxiliary functions which depend on the value of the state over an interval in time, i.e., V(xt), where xt(0) = x(t + 0), —t < 0 < 0.

The conditions for showing an equilibrium point is stable are basically the same as those outlined for the ODE case, above. The main difference comes in showing those conditions are satisfied, which can be more complicated. The Razumikhin approach uses an auxiliary function V(x(t)), but the second condition is relaxed to dV < 0 whenever V(x(t)) > V(x(t + 0)), —t < 0 < 0. Essentially, this requires that V not increase for time intervals longer than the delay.

### 2.4 Phase Models

Many of the analytical tools I have discussed so far are used for studying the stability of equilibrium points and the creation of oscillatory solutions as parameters are varied. These tools are most helpful for predicting the behaviour of systems where the individual neurons do not exhibit oscillatory behaviour when they are uncoupled. For systems which are inherently oscillatory, i.e. systems where the individual neurons exhibit oscillatory behaviour when they are uncoupled, one of the primary tools available is the phase model. The basic idea of this approach is that for a group of oscillating neurons which are weakly coupled, the key variables of importance in understanding how the neurons affect each other are the phases of the oscillators associated with the neurons. Thus a system of k model neurons, each represented by an n-dimensional system of differential equations, can be reduced to a system of k differential equations for the phases of the k oscillators. Typically these equations are in the form where Oi(t) = (9i(t),...,9-i(t),9i+i(t),...,9fc(t)), e = (1,1,..., 1), Q is the network frequency, and e is the strength of the coupling. Since the coupling is weak, e is small, i.e., 0 < e << 1.

The procedure to calculate the phase model for a particular differential equation is described in Hoppensteadt and Izhikevich (1997). In most cases it is not possible to carry out this procedure analytically, however, a numerical implementation is available in the package XPPAUT (Ermentrout, 2005) and described in the book of Ermentrout (2002). The numerical implementation yields a numerical approximation of the functions Hi. A Fourier series representation of these functions can also be calculated.

There are two main results concerning phase models for equations such as (13) which have an explicit time delay in the coupling. The analysis of Ermentrout (1994) and Kopell and Ermentrout (2002) indicates that explicit time delays will produce phase shifts in the corresponding phase models provided that the delay is not a multiple of the oscillation period. Specifically, the models have the form where ^ = tQ mod 2^.

Izhikevich (1998) has refined this analysis. He has shown that Ermentrout's analysis only holds for delays as large as the order of the oscillation period, i.e., t ~ 1/Q. For larger delays, i.e., t ~ 1/(Qe), an explicit delay will occur in the phase model. In this case the phase model will consist of a set of k delay differential equations of the form

For equations with a distributed delay in the coupling, Ermentrout (1994) and Kopell and Ermentrout (2002) have shown that the phase model will be of the form

9i(t) = Qt + e Hi(&i(t - s) - 9i(t)e) ff(s)j ds .

### 2.5 Numerical Tools

There are two basic numerical tools which can aid in the study of delay differential equations such as (13): numerical simulation and numerical bifurcation analysis.

In numerical simulation one attempts to determine an approximate solution of a differential equation given a particular initial state. Note that to solve such a problem for a delay differential equation such as (13), one needs to specify the value of the variable x not just at the start time t = 0, but for the whole interval [—t, 0]. Thus an initial condition for (13) is