I shall consider the situation when two of these neurons are joined with delayed sigmoidal coupling in the following way
•vi(t) = —v3 + (a + 1)v2 — avi — wi + I + c tanh(v2(t — t ) — v) W i(t) = bvi — ywi v2(t) = —v3 + (a + 1)v| — av2 — w2 + I + c tanh(vi(t — t ) — v) w 2(t) = bv2 — yw2
This setup is due to Buric et al. (2005). I will sometimes write (12) in the condensed form x = f (x(t), x(t — t)) , (13)
I will focus on equations with a single discrete delay. The approach is similar for multiple delays, the analysis just becomes more complicated. We will discuss some of the differences that arise for distributed delays in the final section.
There is a very large literature on the effect of time delays on artificial neural networks (ANNs). An example of such a network is the additive (also called Hopfield) neural network with delays. This is usually written in the form n xi(t) = —kiXi(t) + fij(xj(t — Tij) . j=i
I will not attempt to review all the material related to such equations, but will try to highlight those results I feel may have implications for biological neural networks. In particular, networks of the following form n ii(t) = -kiXi(t) + fii(x(t - n)) + E fij(Xj(t - T2) , (14)
j=i have some parallels with biological neural networks, since the uncoupled units may behave as type II oscillators (Campbell et al., 2005).
The main tools for studying the behaviour of delay differential equations are extensions of those for ordinary differential equations which are discussed elsewhere in this volume (Breakspear and Jirsa, 2006). Some familiarity with these tools will be helpful in reading this section.
To improve the flow of the text, I will not give references for all the standard results for delay differential equations that I use. For more information on these, I refer the reader to the fairly accessible books of Kolmanovskii and Nosov (1986) and Stepan (1989) which cover the results of this section or the books of Hale and Verduyn Lunel (1993) and Diekmann et al. (1995) which give complete, although not so accessible, accounts of the theory of delay differential equations.
To begin our discussion, consider the types of solutions which occur most often in neural systems. These are equilibrium solutions (x(t) = x, for some constant X) and periodic solutions (x(t) = x(t + T) for some T > 0). The fundamental questions we would like to answer in order to understand the behaviour of a model with time delays are the following
1. What equilibrium solutions occur in the system?
2. What periodic solutions occur in the system?
3. Are these stable or unstable? That is, do we expect to observe them in experiments and numerical simulations?
4. How do the answers to these questions change as parameters are varied?
Question 1 is easily answered: the equilibrium solutions of a system with time delays are the same as those of the corresponding system with zero delay. Thus for (13) these correspond to constant vectors x such that f (x, x) = 0.
Example. For system (12) the equilibrium points are given by (v1,w1,v2,w2) = (v1,rW1,v2,rW2) where Vj,Wj are constants, found by solving the equations
0 = —V3 + (a + 1)v2 — aVi — W\ + I + c tanh(v2 — v)
0 = —v\ + (a + 1)v\ — av2 — W2 + I + c tanh(v1 — v)) ^ '
It is easy to check that one solution of these equations is (Vi,vj1,V2,vj2) = (v,w,v,w), where v,w are given by (10)-(11). I will focus on this solution in later discussions of this example.
Question 2 is difficult to answer analytically with any completeness. A partial answer can be obtained by determining the bifurcations that occur in the system which lead to the creation of periodic solutions. More detail can be found in subsection 2.2. This question can also be addressed through the use of numerical tools, which are discussed in subsection 2.5.
For equilibrium solutions, question 3 can be addressed via linear stability analysis (see subsection 2.1) and via Lyapunov theory (see subsection 2.3). For periodic solutions this question generally must be answered using numerical tools, as discussed in subsection 2.5.
Answering question 4 is the main goal of bifurcation theory. Analytical methods for studying bifurcations will be discussed in subsection 2.2 and numerical methods in subsection 2.5.
One way to study the stability of an equilibrium solution is through linearization. This is constructed in a similar way as for ordinary differential equations. The linearization of (13) about x is given by x(t) = Ax(t) + Bx(t - t) (16)
where A is a the Jacobian matrix of f (y, z) with respect to y, i.e. the matrix with entries aj = df-, and B is the Jacobian matrix of f(y,z) with respect to z. If the system has multiple delays, then there will be a term in the linearization corresponding to each delay.
It can be shown that, under the right conditions, (16) describes the behaviour of solutions close to x. This will in turn determine the stability of x. To study this behaviour, we assume that there are solutions of (16) of the form x(t) = eAik where A is a complex number and k is an n-vector of complex numbers, to be determined. Substituting this into (16) we obtain
For solutions with k = 0 to exist, we require det[-AI + A + Be-XT]=0 . (18)
If (13) is an n-dimensional system, then (18) can be written in the form n-1
A(A) = An + An-1(Sn-i,0 + Sn-i,ie-XT) + ••• + A ^ 6hje-jXT
where the 5i}j depend on the elements of the matrices A and B.
Equation (19) is called the characteristic equation of the linearization of (13) about x. Any complex number A which satisfies (19) will give rise to a solution of (16) (k is found by solving (17) with the particular value of A substituted in). In practice, we are mostly concerned with the A values for the reasons outlined below.
Example. For our coupled Fitzhugh-Nagumo model (12) the linearization about the equilibrium point (v,w,v,w) is given by (16) where
a -10 0 b -y 0 0 0 0 a -1 0 0 b -y with a = — 3v2 + 2(a + 1) — a, and B =
Note that a depends on all the intrinsic neuron parameters (a,b,j,I), since v is a solution of (10). Putting A, B into (18) shows that the characteristic equation for this example is
Fact: If all the roots of the characteristic equation of the linearization of (13) about xv have negative real part, then xv is asymptotically stable, i.e., all solutions which start sufficiently near to x will tend toward it as t increases.
Fact: If at least one root of the characteristic equation of the linearization of (13) about xv has positive real part, then xv is unstable, i.e., some solutions which start near to xv will tend away from it as t increases.
So we see that to determine the stability of an equilibrium point we need to determine the roots, A of the characteristic (19). These are often called the eigenvalues of the equilibrium point. For ordinary differential equations, the characteristic equation is a polynomial in A and hence there are a finite number of solutions all of which may be calculated or numerically approximated. For delay differential equations, however, the presence of the e-Xr terms means that there are an infinite number of solutions of the characteristic equation. This means we must rely on other methods to determine whether an equilibrium point is stable. Several methods are outlined in the book of Kolmanovskii and Nosov (1986), here we will focus on a particular one which relies on the following result.
Fact: The zeros of A(A) depend continuously on t and the Si,j, and hence on the elements of A and B. Thus as any of these parameters is varied, the number of zeros of A(A) with positive real part can only change if a root passes through the imaginary axis.
The most common way of using this fact in coupled neural systems, is outlined in the following procedure.
1. Set the delay, t, equal to zero. This will change the delay differential equation into an ordinary differential equation with the same equilibrium points as the delay differential equation.
2. Determine the stability of an equilibrium point for the ODE system, i.e. determine the number of eigenvalues with positive real parts.
3. Determine the critical values of the delay, tc < t2 < ■■■ for which the characteristic (19) has eigenvalues with zero real parts. These are the values of the delay where the stability of the equilibrium point may change.
4. Calculate the rate of change of the real part of an eigenvalue with respect to t when t is equal to one of the critical values found in the previous step, i.e., calculate ing, if it is negative, then the number of roots is decreasing. 5. Due to the fact above, the number of roots of the characteristic equation with positive real part will be constant for 0 < t < ti and equal to the number found in step 2. For each subsequent interval, Tk < t < Tk+1, the number can be determined from the number in the previous interval Tk-i < t < Tk, the number of roots with zero real part at Tk and the rate of change calculated in step 4.
Example. Consider our coupled Fitzhugh-Nagumo model (12). We will follow the procedure outlined above.
1. When t = 0 the factors of the characteristic (20) become
2. By analyzing the roots of this equation, it can be shown that if y2 < b def the trivial solution is stable for |c| < j — a = cH, and for c outside this region the equilibrium point has two complex conjugate eigenvalues with positive real part, i.e. it is unstable. (In fact the two points c = ±cH are Hopf bifurcation points for the system with zero delay.)
3. To find the parameter values where the characteristic (20) has eigenvalues with zero real part, we substitute A = iu into (20) and separate into real and imaginary parts. This yields
If dRe(\) > 0, then the number of roots with positive real parts is increas-
Note that we choose the + in the first equation and — in the second for the parameter values for A+ to have a pair of complex conjugate roots and the opposite for A-. Some rearrangement of these equations gives
(bY — a(Y2 + w2))2 + w2(7 2 + w2 — b)2 — c2 (y2 + w2)2 = 0 (21)
Thus, for given values of the parameters a, b, j, I (which determine a) and c one can find w from the first equation and the corresponding t values from the second equation. Alternatively, we can think of these two equations as defining the coupling parameters t and c in terms of the intrinsic neuron parameters and w. Then these equations define curves in the c, t parameter plane. These curves are shown in Fig. 1 for a specific set of intrinsic parameter values. There are multiple curves because tan is a periodic function, i.e., for fixed a,b,Y,w there are multiple values of t that satisfy (22). 4. Taking the appropriate derivatives, we find dX ±\ce
5. Putting together the results of all steps, allows us to fill in the number of eigenvalues with positive real part in each of the subregions of the c, t plane as shown in Fig. 1.
An alternative way to use the procedure outlined above is to set the coupling coefficient (c in (12)) to zero in step 1 and follow the same procedure, varying the coupling coefficient instead of the delay. In systems with multiple delays, the procedure can be followed by setting one of the delays to zero, see e.g. (Campbell et al., 2006, 2005), for examples of this.
To close, we note the work of Olgac and Sipahi (2002, 2005) who have found a way to automate this procedure using a transformation of the characteristic equation.
As noted in the previous subsection, points in parameter space where the characteristic equation has an eigenvalue with zero real part are points where the stability of an equilibrium point may change. These are places where a bifurcation may occur. As discussed elsewhere in this volume (Breakspear and Jirsa, 2006), bifurcations may lead to the creation of other equilibrium points or of a periodic orbit. We refer the reader that chapter for more background on bifurcations.
Was this article helpful?