Bifurcations, stability, and monotonicity properties of a delayed neural network model

PHYSICA Physica D 102 (1997) 349-363

ELSEVIER

Bifurcations, stability, and monotonicity properties of a delayed neural network model Leonard Olien a,b. 1, Jacques Brlair b,c.d,, a Department of Mathematics and Statistics, McGill University, Montrdal, Canada b Centre for Nonlinear Dynamics in Physiology and Medicine, McGill University, Montreal, Canada c D~partement de mathdmatiques et de statistique, Universit~ de Montreal, Montrdal, Canada d Centre de recherches math~matiques, Universit~ de Montreal, Montrdal, Canada

Received 20 May 1996; revised 13 September 1996; accepted 13 September 1996 Communicated by A.V. Holden

Abstract A delay-differential equation modelling an artificial neural network with two neurons is investigated. A linear stability analysis provides parameter values yielding asymptotic stability of the stationary solutions: these can lose stability through either a pitchfork or a Hopf bifurcation, which is shown to be supercritical. At appropriate parameter values, an interaction takes place between the pitchfork and Hopf bifurcations. Conditions are also given for the set of initial conditions that converge to a stable stationary solution to be open and dense in the functional phase space. Analytic results are illustrated with numerical simulations. Keywords: Neural networks; Delay-differentialequations; Hopf bifurcations

1. Introduction The investigation of neural networks has been the subject of much recent activity [ 11 ], to more or less represent the nervous system, or at least to design systems which can perform tasks associated with higher functions of the human central nervous system. One class of models has come to be known as Artificial Neural Networks, since their architecture only metaphorically resemble that of the animal nervous system: pattern recognition problems in particular have been a prime target for these investigations. What came to be known as H o p f i e l d n e t w o r k s have been a prominent tool in the elaboration of these systems. It is of fundamental importance in these studies to determine the convergence of the solutions of a system of differential equations to either one of a number of equilibria, or else to periodic solutions: in the context of pattern * Corresponding author. Dept. Mathrmatiques, Universit6de Montrral, C.P. 6128, Succ. Centre-Ville, Montrral, Qurbec, Canada H3C 3J7; e-mail: [email protected]. I Present address: Division of Medical Genetics, Montreal General Hospital Research Institute, Montrral, Qur., Canada H3A 2T8. 0167-2789/97l$17.00 © 1997 Elsevier Science B.V. All rights reserved PII S01 67-2789(96)00215- 1

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

350

recognition, Content Addressable Memories are nothing more than asymptotically stable stationary solutions with basins of attraction of positive measure. This problem of the determination of boundaries between basins of different attractors has recently emerged as a difficult one, in view of the intricate structure separating adjacent basins. Some of these have been found to have a fractal character [ 17,18] so that structural obstacles to predictability are intrinsically present in some systems. These basin boundaries have recently been investigated mathematically in a model of a neuron with delayed self-connection [20]. In parallel, traditional Lyapunov techniques [5] have been used to establish global stability results. It has been suggested [ 15] and was explicitly recognized in its original formulation [13] that explicit time delays be incorporated in the modelling equations, to take into account the processing time, in each unit, of the information (typically in the form of a voltage) that transits across it. Hopfield's circuit equation [ 13] with no external input but with time delays takes the form

Ciui(t) = - - - u i ( t ) Ri

+

~ j f j ( u j ( t - rj)),

i -- 1. . . . . n.

(1)

j=l

In these equations, ui (t) represents the voltage of unit i which has input capacitance Ci and the connection matrix T = [7~j] has Tij = Rij ( r e s p e c t i v e l y - R i j ) when the noninverting (respectively inverting) output of unit j is connected to the input of neuron i through a resistance Rij : Ri -~- 1/ ~'~j ITij f is thus the parallel resistance at the input of unit i. The transfer function fi (ui) is continuously differentiable, strictly increasing, odd, with limui~ i ~ fi (ui) = 4-1 and f / ( 0 ) > fi' (ui): typically, f (u) = tanh (u), and this is the choice we make (although we expect 'most' calculations to hold for a 'generic' odd function of the same form). Under suitable scaling assumptions [3] (namely, all Ci take a common value C, fi are all the same f = tanh, each unit has equal input resistance) we are lead, for a network with two units, to the system of equations 2

u~(t) = --ui(t) + y ~ a i j f ( u j ( t

- z'j)),

i = 1, 2.

(2)

j=l

This is the system we investigate here. Three common assumptions we do not make at this point are that the connection matrix is symmetric, a12 = a21, that it contains no self-connection, all = a22 = 0, and that units are identical, rl = r2. We eventually somewhat restrict the parameter values, but not uniformly across our results. When initial conditions are specified as continuous functions on the interval [ - max(rj, r2), 0], the basic existence and uniqueness theorems [10] ensure that a unique solution of (2) exists for all positive times. Furthermore, since the function f is odd, it is easy to see that system (2) is preserved under the symmetry (u l, u2) --~ ( - u l, - u 2 ) . We perform a linear stability analysis of Eq. (2) in Section 2; as is usual, stability here means local asymptotic stability. A centre manifold is computed in Section 3, and allows for the determination of the local bifurcations detected by the linear analysis: a degenerate, codimension 2 bifurcation is unfolded into a normal form which leads to an understanding of all possible behaviours when two eigenmodes simultaneously become unstable. In Section 4, we present sufficient conditions for the set of initial functions that converge to an equilibrium to be open and dense in the set of all initial conditions: this result follows from the monotonicity property of the flow generated by system (2).

2. Linear stability analysis We now present the linear stability analysis at the equilibria of (2). Our calculations use the determinant and trace of the connection matrix rather than the eigenvalues of the same matrix, as in [1], or hard to compute parameters

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

351

as in [7]. Also, we show that there is an interaction between a Hopf bifurcation and a pitchfork bifurcation, as in a similar, 3-unit delayed neural network [3,4]. In system (2), the origin (xj, x2) = x = 0 is stationary for all values of the parameters. Other steady states (Xl, X2) satisfy the system

Xl = all tanh(xl) + a12 tanh(x2),

x2 = a21 tanh(xl) + a22 tanh(x2),

(3)

which can be written as 1

x2 = - - [ a 2 2 x j - D tanh(xl)] = f ( x l ) , a12 1

xl = - - [ a l jx2 - D tanh(x2)] = g(x2), a2 I

(4) (5)

where D is the determinant of the network connection matrix. A solution xl of (4) and (5) will thus satisfy xt =

g(f(xl)). It is not difficult to show, using calculus, that nontrivial stationary solutions of system (3) exist provided D < 2T - 1,

(6)

where T denotes one half the trace of the connection matrix. The stability around any equilibrium x = (x~, x~) is determined by the equations .i:l (t) = --Xl (t) + 0"1 lXl (t - rl ) + 0"12x2(t -- l"2), -/2(/) = --x2(t) + 0"21Xl (t - l"1) q- 0"22xz(t -- r2),

(7)

where 0"ij = aij fj (xj). Setting, in the last equations

x(t)-x=eX'(

)c2 cl

(8)

with X a complex number, and cl and c2 two constants, and substituting (8) into (7), we get a nontrivial solution in (cl, c2) if and only if det ()v + 1 - 0",,e-;~rl

--0"21e-)'rl

--0"12e -Xr2 ) )v + 1 -- 0"22e-Xr2 = 0.

This characteristic equation of system (7) determines the stability of the equilibrium solution: the later is stable if and only if all the characteristic roots X, the solutions of ()v + 1) 2 -- ()v -}- 1)(0"11e-'kr' + 0"22e-;vr2) + (0"110"22 -- 0"120"21)e-;v(rt+r2) = 0,

(9)

have negative real part. Finding, in general, all the parameter values for all the roots of Eq. (9) to have negative real parts is a hopeless task: see [16] for analysis of a simpler (!) equation. We consider first a system with both delays equal, rl = r2, and then a network under the common restriction of the absence of self-connection, 0"11 = 0"22 = 0, in Section 2.2.

2.1. Equal delays When rl = r2

=

r,

the characteristic equation (9) becomes

(X + 1)2e 2"~r -- (X + 1)e~r (0"11 + 0'22) + 0"110"22 -- 0"120"21 = 0,

(lO)

L. Olien, J. B~lair/Physica D 102 (1997) 349-363

352

which is a quadratic polynomial in the variable (~ + 1)e xr, and has roots given by (~. + 1)e zr = T 4- v / ~

- D,

(1 1)

where T = ½(u I l + ot22) is half the trace of the connection matrix, and D = a I l ot22 - ot 12ot2 1 is the determinant of the connection matrix. Notice that the expressions on the right-hand side are the two eigenvalues of the connection matrix, a fact used in [I] to investigate networks containing N neurons. We consider in turn the two cases where the connection matrix has real (D < T 2) or complex ( D > T 2) eigenvalues.

2.1.1. Connection matrices with real eigenvalues We write X = p + i09 for a root of the characteristic Eq. (1 0), separate the real and imaginary parts of the ensuing Eqs. (l 1) and obtain ePr [(p + 1)cos(09r) - 09 sin(09r)] = T 4- v / ~

- D,

ePr [09 cos(09r) + (p + 1)sin(ogr)] = 0.

(12) (13)

A change in the stability of the stationary solution can only occur when p = 0, that is cos(09r) - o9 sin(09r) = T 4- V / ~ - D,

(14)

09 cos(09r) + sin(09r) = 0.

(15)

The second equation is equivalent to 09 = - tan(~or). Without loss of generality, we can consider only nonnegative values of 09. If we now write x = 09r, then the roots x of (15) are the points of intersection of the graph of the function y = - t a n ( x ) with the lines y = x / r ( = 09). Thus, co = 0 is a solution for all r , and, if r --~ 0, then (09r) ---~ lzr + k z r , k = 0, 1, 2 . . . Notice that 09 = 0 is always a root of Eq. (15), and that (14) then becomes D=2T-1.

(16)

For values of the parameters D and T on this straight line, the nontrivial stationary states come into existence so that a pitchfork bifurcation occurs at these parameter values. When 09 ~ 0, we can substitute tan(09r) for -09 into (14) to obtain D = 2T sec(09r) - sec2(09r),

(17)

which defines, for a fixed value of r , and for each root 09 of a~ = - tan(09r), a straight line in the plane of the parameters D and T. A m o n g this family of lines, only the one corresponding to values of 09r in the interval (½ Jr, zr) is on the boundary of the stability region, and, when T > sec(09r), this is a line of H o p f bifurcations. At the point of intersection of the lines D = 2 T - 1 and D = 2 T sec (09r) - sec 2 (o~r), there is co-existence of a H o p f bifurcation and a pitchfork bifurcation.

2.1.2. Connection matrices with complex eigenvalues Writing again ~. = p + iw into the characteristic Eq. (10), the separation of Eq. (1 1) into real and imaginary pai'ts now yields ePr[(p + 1)cos(09r) - 09 sin(wr)] = T,

(18)

eP~[wcos(09r) + (p + I) sin(mr)] = 4-v/D - T 2.

(19)

Letting p = 0 leads to

L. Olien, J. Bglair/Physica D 102 (1997) 349-363

353

"i .....................................................

1/3

a

o

113

o

' .......

7" ............

~-- ...........

~ ............

-2

-1

0

~-1

Fig. 1. Stability region (denoted 1) for the null solution in Eq. (10) when r = 1.

cos(oJr) - oJ sin(oJr) = T, co cos(oJr) + sin(~or) = + V ~ -

(20) T2

(21)

so that a change in the sign of o9 does change the sign of the left-hand side of (21). We thus need only consider o9 >_ 0 with + ~ / - D - T 2. Squaring both sides of Eqs. (20) and (21) and adding them give 1 + ~o2 = D ,

(22)

which, upon substituting for ~o in (20), yields T = c o s ( r ~ / D - 1) - v ' - D - 1 sin(rV'-D - 1),

D >_. 1.

Consider now Eqs. (20) and (22) as a pair of parametric equations, parametrized by ~o, for T and D, at a fixed value of r. This curve provides the remaining boundary of the region of stability for Eq. (10), as we shall show. This stability region is determined in the (D, T) plane by the straight lines (16) and (17), and by the parametric curve determined by (20) and (22). It is illustrated in Fig. 1 for the value r = 1. To show indeed that this region actually is the region of stability, we evaluate the derivative of p, the real part of a root of the characteristic equation on each of the three curves forming the boundary of the region, knowing that all characteristic roots have negative real parts inside the region, and only inside the region. Consider first the locus of pitchfork bifurcations given by Eq. (16). This bifurcation corresponds to the connection matrix eigenvalue T + ~ - D, and so if we differentiate, with respect to D, Eqs. (12) and (13), using the + sign in (12), we get, when )~ = p + io9 = 0, Op T

8D

8p +

8D

- 1 --

2~/T 2 - D

so that 8 p / S D is negative. Since D is decreasing as we exit the region along the pitchfork bifurcation curve, the stationary solution loses stability.

L. Olien, J. B#lair/Physica D 102 (1997) 349-363

354

On the Hopf bifurcation curves,we have, for the values D < T 2, the eigenvalue T - x/-T~ - D of the connection m a t r i x . On this bifurcation curve, since m = - tan mr, we have cos(mr) - co sin(mr) • see(mr)

(23)

and so Eqs. (14) and (15) hold using the - sign in (14). Differentiating the latter equations with respect to D, leaves, upon letting p = 0 and rearranging,

Op

cos(mr) - r see(mr)

0 D = 2(T - see(mr))(1 - r 2 secZ(wr))" The Hopfbifurcation occurs when w = - tan(mr) = - sin(mr) sec(mr), with mr ~ (½zr, Jr), so that r sec(mr) = - m r / s i n ( m r ) < - 1 . We know that T - sec(wr) > 0 since T < 1. Thus Op/OD < 0 and, since D decreases upon exiting the stability region along the H o p f bifurcation line, the steady state loses stability when D crosses the line D = 2T sec(mr) - sec2(mr). Finally, on the H o p f bifurcation curve when D > T 2, Eqs. (18) and (19) hold and can be differentiated with respect to D, to give, after some rearrangements, upon substituting p = 0, (cos(mr) -k- r T ) 2

Op sin(mr) + rv/-D - T 2 + sin(mr-----)+ - - ~ ~ J OD

] _

1 24~-Z~"

Here, D increases along the H o p f bifurcation line, and m takes values from zero to wo where o90 is the first positive solution of m = - tan(mr): clearly, wor < Jr, and so, for m c (0, too), sin(mr) > 0. Therefore, sin(mr) +

r ~ / -D- T 2 > 0 and Op/OD > O. 2.2. No self-connection When there is no self-connection in the network, cq 1 =

a'22 = O,

the characteristic equation (9) reduces to

(~. + 1) 2 - ot120~21e-x(rl+r2) = O, which can be factored into (k + 1) = + ~ e

-~(r'+r2/2)

This characteristic equation depends only on the sum of the delays, which we denote 2 r , and allows the reduction to the case analysed in Section 2.1. We have, in the notation of that section, T = 0, and D = -OtlZOt21 so that

D-T2=D. When D < 0, there is a pitchfork bifurcation at D = - 1 . When D > 0, we can use Eqs. (20) and (21) in the case +v/-D, to obtain arcsin (1 / ~/-D) r=

,¢'~- 1

with D > 1 and a stability region which is essentially (use b = - D there) the stability region of a linear scalar delay-differential equation with a single delay [1, Fig. 1].

3. Centre manifold construction We now consider the nature of the bifurcations taking place when the stability of an equilibrium solution of system (2) is lost. Whereas all results of the previous sections concerned linearizafions of system (2) in the form of

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

355

Eqs. (7), we now have to consider the effects of the nonlinear terms of the original equation: these calculations help determine the nature of the bifurcations taking place when the null solution of Eq. (7) loses its stability through an eigenvalue acquiring positive real part, either through the origin or on the imaginary axis. The terminology and notation required for this computation is now introduced, as simply as possible. Consider a system of two functional differential equations written, in standard notation [ 10], as (24)

± = Lxl + g ( x t )

with xt = x ( t + 0 ) , - h < 0 < 0, C = C ( [ - h , 0], g~2), L : C ~

~2 a linear operator, a n d g c c r ( c , ~2), r >_ 1.

L may be expressed in integral form as 0 L05 = / [

(25)

dr/(0)]05(0),

-h where dr/: [ - h , 0] --+ ~2 is a function of bounded variation (a Dirac 3 'function' when a delay differential equation is considered, as is in our case). When Eq. (24) is at such parameter values that the linear part of the equation (26)

x ( t ) = Lxr

has m eigenvalues with zero real parts and all other eigenvalues have negative real parts, then it can be shown that there exists in the state space C an m-dimensional invariant manifold, the centre manifold denoted M r , such that the long term behaviour of solutions to the nonlinear equation is well approximated by the flow on this manifold

[101. For parameter values where the linear equation (26) possesses m eigenvalues with zero real parts, there exists a splitting of the space C = P ~ Q, where P is an m-dimensional subspace spanned by the solutions of (26) corresponding to the m zero real part eigenvalues, and P and Q are invariant under the flow associated with (26). Further, the centre manifold is given by M f = {05 ~ C : 05 = Cbz + h(z, g), z in a neighbourhood of zero in ~m}.

The flow on this centre manifold is Xr = clgz(t) + h ( z ( t ) , g),

(27)

where • (a set of m two-dimensional column vectors represented as a 2 x m matrix) is a basis for P, h c Q, and z satisfies the ordinary differential equation -- Bz + b g ( ~ z

(28)

+h(z,g)).

In (28), B is the m x m matrix of eigenvalues of (26) with zero real part, and b is determined from the solution to the equation adjoint to (26). Indeed, if we let q / ( a set of m two-dimensional row vectors represented by an m × 2 matrix) be the basis for the invariant subspace of the adjoint problem corresponding to P, then b -- tp (0), where tp is normalized by (t/J, qs) : I,

(29)

1 is the m × m identity matrix, and 0 0

05)

(g,(o), 05(0)) - / / @(~ - O)[dr/(O)]05(~) d~ * /

-h 0

(30)

356

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

is the bilinear form associated with (26) and ( , ) represents the usual scalar (dot) product of two vectors. If we let the elements of qJ be linear combinations of those of ~ , i.e. ~ = K~b T (K an m x m matrix of constants), then K = (q[~T, ~ ) - 1 or tp = (~T, q[:~)-l(ib. The problem of describing the long term behaviour of solutions to the two-dimensional system of delaydifferential equations (24) has been reduced (locally) to the problem of describing the behaviour of solutions to the m-dimensional system of ordinary differential equations (28). We now use the construction in this section to first show that the Hopf bifurcation is supercritical, and then unfold the interaction between the pitchfork and Hopf bifurcations. The calculations for this analysis are done using the Maple programs described in detail in [6], and used previously in [2,3]. An analytical, uncomputerized evaluation of the nature of the Hopf bifurcation has also been performed [ 19], but is only sketched here.

3.1. Hopf bifurcations We present first the case of a connection matrix with real eigenvalues. We let rl = r2 = r, D = 0, and T < 0, where D and T are defined in Section 2: since D = 0, the Hopf bifurcation occurs when 2T = sec(ogr) where co = - tan(o9r), o9 6 (lr/2r, z / r ) . Since 2T is the trace of the network connection matrix, this corresponds to a l l + ot22 = sec(ogr). Let us choose oql = - 1 so that a22 = 1 + sec(ogr) = 1 - x, thus defining x. To achieve D = 0 we chooseotl2 = -1,or21 = 1 - t o . In the notation of the beginning of this section, we can compute the basis for the subspace P as

( sin(ogO) • (0) =

ksin(og0)

cos(ogO) ) kcos(og0)

'

where k = [ 0 - + 1)e zr - ~ 1 1 ] / ~ 1 2 . The inner product is gl

0

0

(t/t, cib) = I//1 (0)ci~ 1(0) + 1ff2(0)qb2(0 ) q- Otl l / II11(~ q- "t')~:>l(~) d~ Jr- ot12 / qJl (~ + "()~2 (~) d~ --T

--l"

0 /i +~21 /

0 O2(se + r)q~l (~) d~ + ~ 2 2 /

qJ2(~ + r)q~2(~) d~.

(31)

To compute g' = (q~T, q~)-I as the appropriate basis for the adjoint system, we define C = (q~T, ~ ) , which is a 2 x 2 matrix with entries Cij ~--- (~i, ~j), and compute

-- [(x - 1) 2 + 1][K 2 + 3]

-lv/-~-~2-

1

'

which allows us to evaluate g~ as C - l ~ r. To obtain the ordinary differential system (28) with z =

z2

and

B =

0

'

we consider the power series expansion of tanh, and so

= Bz + bg(dP(-r)z) + O(Izl4), where g contains the cubic terms of f .

(32)

L. Olien, J. Bglair/Physica D 102 (1997) 349-363

357

o

I

I

I

I

-2

-1

0

1

xl

Fig. 2. Phase space representation of an in-phase periodic solution. P a r a m e t e r values are: ~11 = - 1, ot12 = - 2 , ~t21 = - 2 , ot22 = - 3 .

We have b = C-Iq~T(0) and ~i~T(0) = ( 0

0 1

)

x--I

"

The normal form of Eq. (28) can thus be written [9], in polar coordinates, as i" = a r 3 + O(r5),

0 = 09 + O(r2),

(33)

where a is the constant [1 + (x - 1)3][(1 - x 2 ) / 2 -

1]

2 x ( x 2 + 3)

The constant a in the radial equation of system (33) is seen to be negative, since K is positive, and so the Hopf bifurcation is supercritical. This means that a stable limit cycle exists when the equilibrium solution is unstable. Numerically, with r = 1, the smallest strictly positive solution of w = - tan(w) is approximately 2.028758. Since x = - sec(wr), we obtain x = 2.261826 and a = -0.25603. When the condition D = 0 is relaxed, and the connection matrix still has real eigenvalues, an analogous calculation shows that a supercritical Hopf bifurcation takes place: if all = Oil2 = - - 1 , for example, the value of a becomes - 0 . 2 1 2 4 + 0.0303c~22. Since ot22 < - 0 . 2 6 2 8 2 6 along the segment, it follows that a < 0. Considering now a connection matrix with complex eigenvalues, we investigate the loss of stability of an equilibrium along the para .metric curve determined by Eqs. (20) and (22). Computer memory limitations have prevented a consideration of all points of the curve simultaneously, as was done in the case of real eigenvalues of the connection matrix. Five different points were chosen along the length of the curve, and a centre manifold calculation was performed at each point, using in every case ~11 = ~12 = - 1 and r = 1: for each set of parameter values, the constant a was computed to be negative, indicating a supercritical Hopf bifurcation. It is informative to numerically study the periodic solutions born in this Hopf bifurcation, comparing in particular the cases of a connection matrix with real eigenvalues with one in which the eigenvalues are complex. Fig. 2 presents the stable periodic solution in which xl and x2 oscillate in phase: this is the periodic solution created by crossing the

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

358 o e~

u~ c5

tn

I

i

I

t

-2

-1

0

1

xl

Fig. 3. Phase space representation of an out-of-phase periodic solution. Parameter values are: o~11 = 0.3369, ot12 = ot22 = 0.3369.

--

1, or21 = -2.0034,

line D = 2T sec(~or) - sec2(~or). The bold lines indicate the initial function. The parameter values are T = - 2 , D = - 1 which is located at the point ' 2 ' in Fig. 1. In contrast, Fig. 3 shows the stable periodic solution created at the bifurcation corresponding to a connection matrix with complex eigenvalues: here, Xl and x2 are seen to be out of phase. For this illustration, T = 0.336866, D = 2.116851, which corresponds to point '3' in Fig. 1.

3.2. Hopf-pitchfork interaction When the line D = 2T - 1 intersects the curve D ----2T sec(cor) - sec2(ogr) in the (D, T) plane, there coexists a null eigenvalue along with a pair of purely imaginary eigenvalues: a three-dimensional centre manifold can be constructed at these parameter values. The parameter values where this degeneracy occurs are 2T = sec(~or) + 1,

(34)

D=sec(ogr).

(35)

and

Fixing once again r = 1, and otl l ~ - - - O ~ 1 2 ~ - - - 1, then w = 2.028758 and sec(~or) = - 2 . 2 6 1 8 2 6 as in the Section 3.1. Substituting these into (34) gives a22 = - 0 . 2 6 1 8 2 6 and substituting all these in (35) gives a21 = - 2 . 5 2 3 6 5 2 . The system restricted to the centre manifold is, in normal form computed to third order [23],

~=-wy

- (bl(x2d-y2) d - b 2 z 2 ) y + ( a l ( x 2 - b y 2 ) q - a 2 z 2 ) x ,

(36)

p=cox+(bl(x2Wy2)+b2z2)x+(al(x2+y2)q-a2z2)y,

(37)

~ = C l ( X 2 W y 2 ) z q - c 2 z 3,

(38)

L. Olien, J. B~lair/Physica D 102 (1997) 349-363

o

359

~

!_

lal

pitchfork

_

o

~

°

~

zl_ r

Ho Fig. 4. Unfolding of the degenerate flow near the point of intersection of the two different bifurcation curves. The equilibrium points off the z-axis correspond to limit cycles. Secondary bifurcations occur along the lines inside the first quadrant.

where the coefficients al, a2, bl, b2, Cl and c2 depend on the Taylor expansion of the tanh, and the centre manifold itself. This system is best studied in cylindrical polar coordinates with (r, z, 0), since it then takes the form f = a i r 3 + a z r z 2,

(39)

= Cl rZz + czz 3,

(40)

0 =co.

(41)

A computation using the Maple program mentioned before yields the values al = - 0 . 2 2 1 7 8 8 ,

a2 = -1.537082,

bl = - 0 . 0 7 2 8 8 9 ,

b2 = -0.445816,

cl = - 0 . 3 3 0 4 0 9 ,

c2 = -0.429983

(42)

The ensuing system can then be unfolded (see [14]) into f = / z l r - 0.221788r 3 - 1.537082rz 2,

(43)

= / z 2 z - 0.330409r2z - 0.429983z 3,

(44)

0 = 2.028758.

(45)

All possible qualitative behaviours of solutions of the original system (2) which can occur for parameter values near the degenerate bifurcation point must be present in the unfolded system. Fig. 4 shows all possible phase portraits for this three-dimensional system (including the azimuthal coordinate): a rotation about the z-axis must be added to the portraits for a visualization of the full flow. Stationary solutions on the z-axis of the planar system are associated with equilibria of the three-dimensional flow, whereas stationary soltuions off this axis correspond to periodic solutions. As is shown, one of these periodic solutions is stable wherever one exists, except for a single wedge in the first quadrant. A secondary bifurcation is possible, giving rise to a limit cycle not induced by a Hopf bifurcation: this secondary limit cycle is of saddle type and will not directly

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

360

',

'o

;

'2

xl

Fig. 5. M u l t i - s t a b i l i t y d i s p l a y e d in p h a s e - s p a c e : s o l u t i o n s c o n v e r g e to o n e o f t w o e q u i l i b r i u m p o i n t s , in a d d i t i o n to an i n - p h a s e p e r i o d i c s o l u t i o n . P a r a m e t e r v a l u e s are: ot I 1 = - 1, ot12 = 2, ot21 = 3, ot22 = - 1.

affect the observable dynamics. The line labelled ' H o p f ' corresponds to the line, in the parameter plane (D, T), D = 2T sec(wr) - sec2(wr) while the line labelled 'pitchfork' corresponds to D = 2T - 1. Fig. 5 shows orbits originating from different initial conditions (in bold). The parameters correspond to T = - 1, D = - 5 which lies at the point '4' in Fig. 1.

4. Monotonicity results The results of the previous sections were concerned with local properties of the solutions of system (2). The centre manifold construction, although providing information about limit cycles, is fundamentally a local technique, which traces where those periodic solutions are born in parameter space. In contrast, in this section, we apply a totally different set of techniques, topological, which yield global information on a set of solutions. For this, we go back to system (2) to which we apply the notation introduced in Section 3, allowing us to write ft" (~0) = --~Oi (0) + a i l tanh(~oi ( - r l ) ) + ai2 tanh(~o2(-r2)),

i : 1, 2.

For a system of differential equations inducing a strongly monotone flow on its phase space, the set of initial conditions that converge to a stable equilibrium is known [ 12] to be an open, dense subset of the phase space. In particular, any numerically computed trajectory is expected to converge to an equilibrium. A sufficient condition for a system to generate an eventually strongly monotone flow [22] is that it be cooperative and irreducible. We present a proof, along the arguments of [22], that for some parameter values of the connection matrix entries otij, system (2) is strongly ordered. We first recall the necessary concepts. In the space of continuous functions C r = C ( [ - r l , 0], R) x C ( [ - r 2 , 0], R ) , we consider the partial ordering induced by the subset C~_ = {~o(0) = (x(O), y(O)) ~ C r Ix(0) ___ 0 V 0 ~ [ - r l , 0], y(O) > OVO ~ [ - r 2 , 0]}: we write tp < ~p iff ~ - tp ~ C~_, ~o < ~p iff ~o _< ~p and ~o # ~p and we write ~0 << ~p iff lp - tp ~ Int C~_ and ~0 # ~p. We show that system (2) is eventually strongly monotone, and hence monotone, provided a21al2 > 0 and all, a22 >_ O: the set of initial conditions converging to a stable equilibrium is then open and dense in C r. In

361

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

the notation of Section 2, T = l ( a l l -I- a12) >_ 0 and D = alia22 - a12a21 _< T 2 under these conditions: the region T > 0 and D < T 2 does not intersect the regions of the parameter plane (D, T) where periodic solutions were shown to exist in the last section. They provide, in a complementary fashion, a set in parameter space where, essentially, the only stable sets are equilibria. System (2) induces a flow 4~ in the phase space C r in the usual way: for xo c X define 4~t(xo) = x(t) where t ~ x(t) is the solution to the system of differential equations such that x(0) = xo. The flow ~b is said to be (strongly) monotone if the map q~t is such for each t > 0; similarly, ¢ is eventually strongly monotone if there exists to < o~ such that 4~t is a strongly monotone map for all t _> to. Our precise result is:

Theorem 1. Assume a l l , a22 >_ 0 and ot12, or21 > 0. Suppose ~o, ~p E C r and ~0 < ~ . Let us denote by x(t, ~o) the value of the solution, at time t, of system (2) generated by the initial condition (function) qS; then x(t, q~) << x(t, ~p) for t > rl + r2. The result follows (see proof below) from the following two lemmas.

Lemma 2. For any two (vector) functions 4~ and ~p in C r satisfying ~0 < ~ , ~oi(0) = 1//i (0) for some i I1 or 2), then fi(~o) <_ fi(~P) (for the same i) Proof We have fi (lp) - fi (qg) = all tanh(apl ( - r l ) ) + ai2 t a n h ( ~ 2 ( - z 2 ) ) - ail tanh(~01 ( - - r l ) ) -- ai2 tanh(~o2(-r2)) >_ 0 since a l l , a22 > 0 and a l2, a21 > 0 and the hypotheses on the order of ~p and ~p.

[]

Lemma 3. Suppose, in the notation of the previous lemma, that ~o < 7t. Then, in the notation introduced before the statement of the theorem, x(t, ~o) < x(t, ap), V t > O. Proof Let fE(~P) = f(~P) + e(1, 1) T, and ~pe = ~ + e(1, 1) T. We will show that x(t, qg, f ) < x(t, ~e, fe), then let e ~

0 and invoke the theorem on continuous dependence on parameters. Suppose the proposition is false. Then

there exist e > 0 and t' > 0 such that x(t, ~p, f ) < x(t, ~ , fe) on I - r , t') and xi(t', ~0, f ) = xi(t', ~ , fe). This implies that 2i(t', ~o, f ) >_ xi(t', ~re, fe). But 2 i ( t ' , ~e, fE) = f i ( x t ' ( ~ , fe)) + e > fi(xt'(Oe, f~)) > fi(xt'(g o, f ) ) = 2i(t', ~p, f ) which is a contradiction.

[]

Proof of theorem 1 We have x(t, ~o) < x(t, ~ ) for all t _> 0. We can write 1

X(I,

x ( t , qo)

] d~ox(t, sTt + (1 - s)~p)(~ - cp) ds q/

0 and it suffices to show that the integrand is strictly positive for t > rl + r2. Let se c C r, ~ ~ C~_ - {0} and define y(t, ~) = d¢x(t, ~)~. Then y satisfies the linear functional differential equation

f(t) = d f (xt(~))yt with initial condition Y0 ~ f l .

362

L. Olien, J. Bdlair/Physica D 102 (1997) 349-363

Let L(t, 4) = d f ( x t ( 4 ) ) then, for ( e C r, 0!11(1 (__Z,l) - - ( l (0) + cosh2(xl (t - r l ) ) -+- cosh2(xz(t - "r2)) L(t, 4)( = 0/21 (l ( - ' t l ) 0~22(2(-r2) ~k--(2(0) -k- cosh2(xl (t - r l ) ) -k- cosh2(xz(t _ r2))

=

-(1(0) - ( 2 (0)

+ L(t,4)(.

Assume that y l (t') > 0. Now fl > 0, and it follows from an invariance result [21 ] that Yt > 0, so that L l (t, ~)yt > 0, for any t > 0, and hence y l ( t ) = - y l ( t ) + Ll(t, 4)Yt >_ - - y l ( t ) . Thus, for t > t', Yl (t ) >_ Yl ( t')e -(t-t') > 0. A similar calculation holds for Y2 ( t'), and so we obtain that if yi (t') > 0 for some t' > 0 then Yi (t) > 0 for any t > f . We now show that Yt > 0 for t > rl + r2. First assume/3(0) = 0. Since/3 E C~_ - {0}, there exists i and t ' e [ - r i , 0) such that/3i(t I) > 0. Recall that Yt > 0 for t > 0. Without loss of generality, assume/31(t') > 0. If

y2(t' + r l ) = 0 then, since yi(t) =/3i(t) for t E [ - r i , 0), i = 1,2, we have t~21/31(t') ~22Y2(t' + rl -- r2) ~2(t' + r l ) -- cosh2(xl(t,) ) + cosh2(x2(t r + rl - r2)) > 0 and this is a contradiction because it would imply that every left neighbourhood of t' + rl would contain a point t* with y(t*) < 0. So yz(t I + r l ) > 0 and hence yz(t) > 0 whenever t > t t + r l . Now suppose that Yl (t r + rj + r2) = 0. Then OtllYl(t' + r2) OtlZY2(t' + r l ) + >0, Yl (t' -k- rl + r2) -----coshz(xl (t' q-- rz)) coshZ(xz(F -k- r l ) ) which again is a contradiction and thus establishes Yl (t) > 0 whenever t > t' + rl + r2. Now if one of/31 = (0) or/32 : (0), then the first part of the argument can be repeated with t' = 0. Finally, if /3(0) -¢ 0, then the result follows directly from the manipulation above.

[]

Corollary 4. If a l l , a22 > 0 and a12, a21 < 0 in system (2), then the set of initial conditions that converge to a stable equlibrium is open and dense in C r. Proof Let x3 = - X l , a23 = - a 2 1 , a32 = - a 1 3 , a33 = a l l and r3 ---- r l . By the symmetry of system (2), the new system in the variables x2, x3 meets the conditions of the main theorem.

5. Discussion We have studied the stability in a model system for a neural network, when time delays are incorporated to take into account the transmission of the signal inside the units. Bifurcations were detected, and their nature was elucidated by the construction of a centre manifold, a reduction to normal form, and an unfolding of bifurcations of co-dimensions one and two. Regions of parameter space where multistability is possible were identified, as well as co-existing limit cycles. This work complements that of Gopalsamy and He [7] and Gopalsamy and Leung [8]. In the first paper, the authors consider a system of dimension n for which Lyapunov functionals are constructed to show that an equilibrium is globally asymptotically stable, for all values of the time delays. This result is obtained by restricting somewhat

L. Olien, J. B~lair/Physica D 102 (1997) 349-363

363

both the values o f the entries o f the c o n n e c t i o n matrix and the m a x i m u m slope of the transfer function f . In the second paper, a t w o - d i m e n s i o n a l system without self-connection is analysed, for its linear stability and also for the presence of H o p f bifurcation. The nature of the latter is determined by a perturbation expansion. As in previous w o r k [3], we have no particular need for the s y m m e t r y o f the c o n n e c t i o n matrix. Such an assumption w o u l d m e r e l y simplify our analysis (the eigenvalues w o u l d always be real, for example). The most interesting aspect o f our work concerns the existence, and stability, o f invariant sets m o r e c o m p l i c a t e d than equilibrium points, namely limit cycles. Further, one limit cycle does not arise directly from a H o p f bifurcation, but rather from the interaction with the transcritical bifurcation. For the more interesting bifurcations to o c c u r in our system, the networks need to possess self-connections. This assumption differs s o m e w h a t from a c o m m o n hypothesis in most neural network architecture [ 1 1].

Acknowledgements The authors a c k n o w l e d g e financial support to JB f r o m the Natural Sciences and E n g i n e e r i n g Research Council (Canada) and the Fonds pour la F o r m a t i o n de Chercheurs et l ' A i d e ~ la R e c h e r c h e (Qu6bec).

References Ill J. B61air, Stability in a model of a delayed neural network, J. Dynamics and Differential Equations 5 (1993) 607-623. 121 J. B61air and S. Campbell, Stability and bifurcations of equilibria in a multiple-delayed differential equation, SIAM J. Appl. Math. 54 (1994) 1402-1423. 131 J. B61air, S. Campbell and P. van den Driessche, Frustation, stability and delay-induced oscillations in a neural network model, SIAM J. Appl. Math. 56 (1996) 245-255. 141 J. B61air and S. Dufour, Stability in a three-dimensional system of delay-differential equations, Canad. Appl. Math. Quarter. 4 (1996) 135-156. 151 T.A. Burton, Neural networks with memory, J. Appl. Math. Stochastic Anal. 4 (1991) 313-332. 161 S. Campbell and J. B61air, Analytical and symbolically assisted investigation of Hopf bifurcations in delay-differential equations, Canad. Appl. Math. Quarter. 3 (1995) 137-154. 171 K. Gopalsamy and X.Z. He, Stability in asymmetric Hopfield nets with transmission delays, Physica D 76 (1994) 344-358. 181 K. Gopalsamy and I. Leung, Delay-induced periodicity in a neural netlet of excitation and inhibition, Physica D 89 (1996) 395-426. [91 J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer, New York, 1983). 110] J. Hale and S.V. Lunel, Introduction to Functional Differential Equations (Springer, New York, 1993). [ 111 J. Hertz, A. Krogh and R. Palmer, Introduction to the Theory of Neural Computation (Addison-Wesley, Redwood City, CA, 1991 ). [I 21 M. Hirsch, Stability and convergence in strongly monotone dynamical systems, J. Reine Angew. Math. 383 (1988) 1-53. [ 13 ] J.J. Hopfield, Neurons with graded resonse have collective computational properties like those of two-state neurons, Proc. Nat. Acad. Sci. 81 (1984) 3088-3092. [ 14l W.E Langford, A review of interactions of Hopf and steady-sate bifurcations, in: Nonlinear Dynamics and Turbulence, eds. G.I. Barenblan, G. Iooss and D.D. Joseph (Pitnam, Boston, 1983) pp. 215-237. 115] C.M. Marcus and R.M. Westervelt, Stability of analog neural networks with delay, Phys. Rev. A 39 (1989) 347-359. [16] J.M. Mahaffy, K.M. Joiner and P.J. Zak, A geometric analysis of stability regions for a linear differential equation with two delays, Internat. J. Bifurcations and Chaos 5 (1995) 779-796. [ 171 S.W. McDonald, C. Grebogi, E. Ott and J. A. Yorke, Fractal basin boundaries, Physica D 17 (1985) 125-153. 1181 H.E. Nusse, J.A. Yorke, Wada basin boundaries and basin cells, Physica D 90 (1996) 242-261. [191 L. Olien, Analysis of a delay differential equation model of a neural network, M.Sc. Thesis (Mathematics), McGill University, Montr6al. [20] K. Pakdaman, C.P. Malta, C. Grotta-Ragazzo and J.-F. Vibert, Effect of delay on the boundary of the basin of attraction in a self-excited single graded-response neuron, Neural Computation (1997), in press. [21 ] G. Seifert, Positively invariant closed sets for systems of delay differential equations, J. Differential Equations 22 (1976) 292-304. [22] H. Smith, Monotone semiflows generated by functional differential equations, J. Differential Equations 66 (1987) 420-442. 1231 S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos (Springer, New York, 1990).