1
Multistability in Coupled Fitzhugh-Nagumo Oscillators
Sue Ann Campbellab*t and Michael Waited
aDepartment of Applied Mathematics, University of Waterloo Waterloo, Ontario N2L 3G1.
bCentre for Nonlinear Dynamics in Physiology and Medicine, McGill University Montreal Quebec.
We consider a pair of neurons modelled by Fitzhugh-Nagumo equations with electrical coupling. When the neurons are identical, we show how the symmetry of the system leads to the coexistence multiple, stable periodic orbits. As the coupling between the neurons is strengthened, these periodic orbits can undergo various bifurcations, leading to the coexistence of multiple, stable chaotic attractors. We show that this behaviour persists when the neurons are close to, put no longer exactly, identical.
1. Introduction
The Fitzhugh-Nagumo [1,2] equations are an set of simple equations which exhibit the qualitative behaviour observed in neurons, namely: quiescence, excitability and periodic behaviour. The form we will use here
x3 1
x = c(y +x - — + z(t)), y---(x-a + by), (1)
■i c
was introduced in [1] by modifying the equations of the van der Pol relaxation oscillator [3]. Although the variables have no exact physiological interpretation, for appropriate parameter values (see section two), the qualitative behaviour of x is similar to that of the voltage variable in the Hodgkin Huxley equations [4] and that of y to the "recovery" or gating variables. The function z(t) represents the forcing of the cell by an external stimulus. We will be interested in the case when there is no external forcing, z(t) = 0, so we may focus on effects of coupling two such neurons together. In section two we will review the possible behaviour which may occur in (1) as the parameters a, b, c are varied. This will allow us to choose reasonable range of values for our coupled neuron studies.
Our study is designed to emulate two neurons linked with electrical coupling, i.e. coupling via the flow of ions through the gap junctions between neurons. We thus follow [5] in choosing the form of this coupling to be a constant times the difference in the voltage of the two cells. As in this work (and differing from [6]) we couple the neurons only through their voltage equation. This leads to a set of four coupled, nonlinear ordinary differential equations:
x3 1
«i   =   ci(t/i + xx - -3-) + 7i2(«i - x2),   ýi   =--(xx-ax + bxyx),
O Cl
X3 1
«2    =    C2(t/2 + X2 - -£■) + 721 («2 - «l),    2/2     =--(«2 - 02 + &2Í/2)-
(2)
C2
From a biological standpoint, if we consider the two neurons to be in a similar region of the brain, it is likely that the parameter values will be similar but not identical. Thus our focus will be on the case where the parameter values are such that both neurons will be capable of exhibiting the same qualitative behaviour.
*To whom correspondence should be addressed. Email: sacampbell@uwaterloo.ca
t Supported by the Natural Sciences and Engineering Research Council of Canada and the MITACS Network of Centres of Excellence.
Í Current Address: Department of Meteorology, McGill University, Montreal Quebec.
2
To give us a basis from which to start, we consider, in section three, the behaviour when two neurons with identical parameters are coupled together (i.e. (2) with a\ — a2, b\ — b2, c\ — c2 and 712 — 721)-We will focus on parameter values such that the neurons may display either a single stable steady state or periodic behaviour. This will lead us to section 4, where we study the case of near-identical neurons, i.e neurons where the parameters differ by a small amount. In section 5 we discuss our results in relation to some other similar studies.
2. Single Neuron
Consider equation 1 with no external forcing
x — c(y + x ■
1
V----a + by).
Equilibrium solutions, (x(t), y(t)) — (x, y), must satisfy
x3 + 3 ( -
1 ) x - 3- = 0,
y
(3)
(4)
for 6^0 and [x, y) — (a, a3/3 —a) when 6—0. We will restrict our analysis to b > 0 since the expressions for the equilibrium points are well behaved in this interval, and since this will be the relevant parameter range for the studies of subsequent sections.
One, two or three solutions to (4) can exist, depending on whether the quantity D — a2 +4(1 — b)3/9b is positive, zero or negative, respectively. Verifying the standard conditions it can be shown that a saddle node bifurcation occurs along the set of parameter values where D — 0.
The Jacobian of linearization of (3) about an equilibrium is
c
which has eigenvalues
(5)
c(l — X ) — — i c
y/(c(l-^)_^2 + 4(6(l-^)_1)
Therefore, a Hopf bifurcation of the equilibrium point (x, y) can occur when b2 < c2 and
x' -l + — = 0.
Using (4) shows that this will occur when b2 < c2 and
b
b
X~~2
b (      b ,
3
a = 0.
(6)
(7)
We will focus on parameter values where only one equilibrium point exists, but this point may undergo a Hopf bifurcation. This is the parameter range chosen for biological reasons by Fitzhugh [1]. Clearly this corresponds 0 < b < min(l, |c|). For any fixed b in this range there can be up to two Hopf bifurcations of the equilibrium as a is varied.
3. Coupled, Identical Neurons
We now turn to the situation where two iden case is (2) with cij — a, bj — b, Cj — c; j — 1, 2. 712 — 721 — 7, yielding the equations
x3
xi - c(j/i + xx - -^-) +7(^1 - x2), yi = «2   =   c(j/2 + x2 -      +f{x2 -        y2 -
ical neurons are coupled together. The model for this We further assume that the coupling is symmetric, i.e.
— (xx -a + byx),
I (8)
— (x2 -a + by2).
3
Note that a consequence of the form of the coupling and the assumption that the neurons are identical is the invariance of the equations under the transformation [x±, y±, X2, 2/2) (*25 2/25 J/i)- This symmetry can also be observed in the existence of an invariant subspace for the equations I — {(xi, j/i, X2,2/2)1*1 — X2,yi — 2/2}- We will refer to solutions which lie in I as symmetric and those which do not lie in I as non-symmetric solutions. Physically, the invariance of this subspace means that if the two neurons start with identical initial conditions then their subsequent behaviour will also be identical; such solutions are sometime referred to as in-phase.
To facilitate our study of these equations we will perform a simple, linear change of variables
%2)
■2/2)
X2 = \{%i + «2)
Y2 = i (2/1 + 2/2)
leading to the new equations Xl   = c(Y1+X1-^--
X1Xi)+2yX1,
X2 = c{Y2+X2-^--X(X2), In the new variables the symmetry is (X1,Y1,X2,Y2) (-X1,-Y1,X2,Y2)
Y2
— (Xí + bYí), c
--{X2-a + bY2). c
(9)
(10)
(11)
and the invariant subspace (and hence the symmetric solutions) have the simple form I — {(Xi, Y±, X2, Y2) \X± — Y± — 0}. A further symmetry is also evident
(x1,y1,x2,y2;«) (Xi,Yi,-X2,-y2;-a).
(12)
Once again, we begin our study by determining the equilibrium points of the equations. Symmetric equilibrium points are of the form
{XuYuX2,Y2) = {0,0,x, y)
where x and y satisfy (4). Non-symmetric equilibrium points are given by
(Xi, Yi, X2, Y2) = (±X1,±Y1, X2, %)
where
Yi -
1
Yn =
a-X2
and X2 satisfies
1 "
1
l-b
+ 9l c
3a
X2 + 86-°
(13)
(14)
(15)
(16)
The Jacobian of the linearization of (10) about the symmetric equilibrium is
3(x) =
c(l - x2) +27 _ 1
c
0 0
0 0
c(l-x2)
(17)
Clearly, the four eigenvalues of this matrix come in two pairs.
The symmetric eigenvalues have eigenvectors lying in the invariant subspace I and satisfy A2 — /JjA + aj — 0, with Px — c(l — x1) — b/c and ax — 1 — b{l — x1). Bifurcations associated with these eigenvalues will result in new solutions which remain in I. We will refer to these as symmetry preserving bifurcations. These bifurcations are just the bifurcations of the Fitzhugh-Nagumo equations studied in the section two. Recall that the saddle node bifurcation occurs on
a2 + — (1-
b3) - 0
4
while a Hopf bifurcation can occur on 6
1 -
+ 1-6
- a = 0,   if bl < c\
(18)
The non-symmetric eigenvalues have eigenvectors which lie in the complement of I and satisfy A2 — AvA + ctM — 0 with — c(l — x1) — b/c + 2j and a_v — 1 — 6(1 — x1) — 26/07. Bifurcations associated with these eigenvalues will result in new solutions which no longer lie in I. We will refer to these as symmetry breaking bifurcations. Proceeding in a similar manner to the last section, we find that there can be a Hopf bifurcation when fitf — 0 if atf > 0, that is when
^_1 + ^_^7=0 if62<c2_
cl c
Using the expression for the equilibrium points (4) this becomes
cz c
cz c
+ 1-6
0,
if 62 <c2.
Further, there can be a steady state bifurcation when ctjj — 0, which upon using (4) becomes 1
1 2
- + -7 6 c
1
b + h
+ 1-6
a = 0.
(19)
(20)
(21)
Consideration of the form of the equilibria (14) (which is a consequence of the symmetry (11) of the equations), we see that this must be a pitchfork bifurcation.
Note that when 7 — 0, both (18) and (20) will be satisfied simultaneously. At such parameter values, the
equilibrium point will have a repeated pair of pure imaginary eigenvalues: A — ±«y^ 1 — Consideration of the Jacobian of the linearization at such parameter values,
M 0
0 M
with
M
r - c 1	0 =	' 0	0 "
c          c .		0	0 _
(22)
shows that these are points of D2 symmetric Hopf bifurcation. As shown in [16], such points can be the source of more complicated dynamics resulting from secondary bifurcations.
Now consider the non-symmetric equilibria. Clearly, there can be one, two or three solutions of (16). The transition from one to three corresponds to a saddle node bifurcation of the non-symmetric equilibrium points. This will occur when the discriminant of (16) equals zero, which occurs when the
+
0.
(23)
The Jacobian of the linearization of (10) about the non-symmetric equilibrium points is " c(l - X\ - X2) + 27 c
J(xl5x2)
-2cXiX2 0
-2cXiX2 0
c(L
and the corresponding characteristic equation is
6 c
A2 +
c    2c -,       \ ,    26 -, 26 ; + y*2 + 27   A+yX2 + -7
A2 +
+
2c -.
-AbcXt
1 1
27
1----Xi + -!-
b    3   1 c
A +
A+^
c
26
0.
To determine parameter values where a Hopf bifurcation may occur, we let A — iu in this equation, separate into real and imaginary parts and eliminate ui. This gives a constraint on the parameters for the existence of a Hopf bifurcation: ^4(7; a, 6, c)X\ + ^3(7; a, 6, c)X2 + ^5(7; a, 6, c) — 0, where Pj(~f, a, 6, c)
5
0.14
Figure 1. (a) Bifurcation sets in 7, a plane for coupled identical oscillators (10). Here b — 0.4 and c — 2, but the diagram would be similar for any b, c satisfying b < min(l, |c|). (b) Bifurcation diagram of Xi vs 7 for the same parameter values and a > Hs fixed (a — 0.7). Solid (dashed) lines correspond to stable (unstable) equilibrium points, closed (open) circles correspond to the maximum amplitude of stable (unstable) periodic orbits.
is a j order polynomial in 7 with coefficients which are functions of a, 6, c. Using the equilibrium point equation, we can then derive an implicit Hopf bifurcation set involving 7, a, b, c. Due to its length, we do not reproduce the expression here.
We illustrate our results by fixing b and c and drawing the curves corresponding to the bifurcation sets described above in the 7, a plane (Fig. 1). We have chosen the case where 0 < b < 1 < c. Thus in the absence of coupling (7 — 0), each neuron can exhibit either steady state or periodic behaviour. The labelling of curves is as follows Hs - symmetry preserving Hopf bifurcation, Hb - symmetry breaking Hopf bifurcation, H„ - Hopf bifurcation of non-symmetric equilibrium point, P - pitchfork bifurcation, SN„ - saddle node bifurcation of non-symmetric equilibrium point. Recall from the analysis of section two that the saddle node of the symmetric equilibrium point cannot occur for these parameter values.
Numerical investigations near the symmetric Hopf bifurcation points (the intersections of the lines Hs with the curve Hb in Fig. 1(a)) reveal that secondary bifurcations arise from these points. For example, the bifurcation diagram of Fig. 1(b) shows period doubling and pitchfork bifurcation of limit cycles. These bifurcations lead to the coexistence of multiple stable limit cycles (Fig. 2(a)) and of a stable limit cycle with more complicated attractors (Fig. 2(b)).
4. Coupled, Near-Identical Neurons
It is clear from the last section that the symmetric Hopf bifurcations which resulted from using identical models for the two neurons can be a source of multistability. It should be expected that when the neurons are no longer identical this symmetric Hopf bifurcation point will no longer exist. A little reflection reveals two possible situations in the non-symmetric case: (a) the symmetric Hopf bifurcation point is replaced by a non-symmetric Hopf-Hopf interaction point (b) there is no Hopf-Hopf interaction point at all. We show below that the former occurs if only the coupling parameters are different for the two neurons, but the latter occurs if just one of the parameters differs between the equations modelling the two neurons.
The simplest case occurs when all the parameters of the neurons are the same, only the coupling parameters differ. The model for this case is (2) with cij — a, bj — b, Cj — c; j — 1, 2. For comparison
6
Figure 2. Numerical simulations of (10) with a — 0.7, b — 0.4 and c — 2 for three initial conditions resulting in three different attractors. (a) Attractors are all limit cycles for 7 — 0.086. (b) Two attractors have undergone period doubling bifurcations by 7 — 0.087.
with the symmetric case we also set 712 — 7 and 721 — 7 + e, yielding
«1   =   c(j/i + xx - ^-) + f(x1 - x2)
„3
2/1
X 9
x2
c(j/2 + «2 - y) + (7 + <0(«2 - «1), £2
- a + 6j/i),
c
--(x2 -a + by2).
c
(24)
Applying the transformation (9), these become Xi   -   c(y1+X1-^-X1X22) + (27 + e)X1,   Yi   - +feY^,
X33
X2   - c{Y2+X2-^--X21X2)-eX1,
c
Y2   =   --(X2-a + bY2). c
(25)
The analysis of this case is very similar to the symmetric case.  In particular, the subspace I is still invariant and the symmetric fixed points are given by (13) and (4). The Jacobian of the linearization about these fixed points is
J(*) =
c(l - x2) + 2f + e c 1 b
—e 0
0 c(l - x2 0 -h
(26)
Similarly to the symmetric case the eigenvalues split into those with eigenvectors in I and those with eigenvectors not strictly in I. Proceeding as in section three, we find that the symmetry preserving Hopf bifurcation can occur on
±W1-
b (      b , 3 1-^1+1-6
- a = 0,
if 62 < c2.
and the symmetry breaking Hopf bifurcation on b
±W1_4 + ?7±£
1-^ +
27 + e
+ 1-6
- a = 0,
if 62 <c2.
(27)
(28)
Thus when 7 — — e/2, equations (27) and (28) may be satisfied simultaneously. By contrast with the symmetric case, the Jacobian of the linearization at such points in this case is
M 0 E M
with M
r - c 1	,     E =	' -e   0 '	.     0 =	' 0	0 "
c _ 1 _b c          c .		0 0		0	0 _
(29)
7
W4
-4
	
	
	
	
0 5	
(b)
1.5
7
-1-1-1-1— ■   y^' '    ' '	if ° cb .          .      . •        mo   <bo °
- V • 1           1           1 1	.....
0.1
0.3
Figure 3. (a) Bifurcation sets in 7, a plane for identical oscillators with non-identical coupling (25). Here b — 0.4, c — 2 and e — —0.2, but the diagram would be similar for any b, c satisfying b < min(l, |c|) and ^gjjL 0. (b) Bifurcation diagram of X\ vs 7 for the saijt^ parameter values and a > Hs, fixed (a — 0.7).
Figure 4. Numerical simulations of (25) with a — 0.7, b — 0.4, c — 2 and e — limit cycle coexist at 7 — 0.1337. (b) Two limit cycles coexist at 7 — .24
-0.2. (a) A torus and a
Hence the points correspond to one-to-one resonant Hopf-Hopf interactions. Resonant Hopf bifurcations have been studied in [7], where it is shown that such secondary bifurcations as saddle node bifurcation of limit cycles and Neimark-Sacker bifurcations can result.
The other bifurcations, including those of the non-symmetric fixed points, can be analyzed using a similar approach as was used in section three. However, since we are primarily interested in the presence of multistability in these equations, we don't present this analysis here. Instead we show a numerically generated bifurcation set for (24) in the 7, a plane for b — 0.4, c — 2 and e — 0.2 in Fig. 3(a). The labelling of curves is as in Fig. 1(a). Note that the intersection points of the symmetry preserving and symmetry breaking Hopf bifurcations have been shifted to the right as predicted by our analysis.
Fixing a at a value slightly above Hs and varying 7 close to the intersection point reveals secondary bifurcations (Fig. 3(b)), which lead to multistability (Fig. 4).
To study what happens when the parameters in the models for the two neurons are different, we consider the case where all the parameters in the equations for the neurons are the same except for the Cj. We take this case in particular as maintaining the equality of the a's and 6's means that this model
8
will still admit symmetric equilibria, which renders the analysis simpler. It also means that this is the "least harmful" way of breaking the symmetry of these equations. The resulting model is
X\     =     Cl(j/l + Kl
x2
-    C2(j/2 + X2 - —) + 7(«2
X2)
X\)
2/1
2/2
(x1-a + 6j/i)
(30)
---(x2-a + by2)
c2
One could also write equations for the transformed variables (Xi, Yi, X2, Y2) defined by (9), however these are considerably more complicated than (30). For comparison with section three, we will use the transformed variables to represent solutions graphically.
A brief analysis of (30) reveals that although the subspace I is no longer invariant, the symmetric equilibria still exist and are identical to those of (8), i.e. they are given by [x±, y±, x2, y2) — (x, y, x, y), with x, y solutions of (4). The Jacobian of the linearization of (30) about these equilibria is
3(x)
ci(l
-x2) +7
-7 0
ci
_b_
0 0
-7 0
c2(l - x2) +7 _ j_
C2
C2 _b_
C2
(31)
where x is a solution of (4). The characteristic equation for this system is then
a4 +
(ci+c2)(k2-1) +
b (ci + c2) cic2
27
a3 +
2 +
(clC2(x2 - 1) + 6)2
+b2±±Sl(x2
C\C2
(ci +c2)(x2 - 1) +
(ci +c2)(x2 Cl +c2
1)
26 (ci + c2)
7
+ (b{a
1) + 1
cic2 ft (ci + c2
ClC2
ClC2
_ (6(2+^+^) (^_i) + 2(—+ 1
cic2 y Vcic2
ClC2
a2 + M ' b2
x2
1) + 1) X
7 a
(32)
(6(
x2
1) + 1) 7 = 0
Clearly, when 7 — 0, the characteristic polynomial factors into two quadratics, yielding two characteristic equations a2 + (b/cj + Cj(x2 — 1)) a + 1 + b(x2 — 1) — 0, one for each neuron. For fixed b and ci ^ c2, the Hopf bifurcations of the two uncoupled neurons now occur at two distinct values of a
±Wi-
1 -
+ 1-6
def
OHj,
(33)
When 7 is nonzero, one might expect that there will still be two sets of parameter values along which a Hopf bifurcation occurs, corresponding to the two sets above. To find the expressions for the sets of parameters along which a Hopf bifurcation can occur, we put a — iu into the characteristic equation (32) and separate the resulting equation into real and imaginary parts. Isolating 7 yields
7 - (ci +c2)(c1c2(x2 where ui satisfies
1)+*)
b{x2 - 1) + 1 - J1
b(x2 - l)(ci + c2)2 + 262 + 2cic2(l - ui2
co6 + (c2 + 4)
+2b(x2 - 1) +
62
0      '-'2      o^„2  , „2\
- 3 +
ci +ci
2c2c2
{2b2 - 3{ci + c22)) + 3
r2r2 clc2
(ci + 4)-i (x'-i
J2 + (b(x2 - 1) + 1)
b2(cj + c22) '^c\c2
-1
(34)
(35)
0.
Thus solving (4) for x and (35) for ui and substituting in (34) yields the desired expression for 7 as a function of a, b and the Cj. This expression is quite lengthy, so we don't reproduce it here, instead we show a graph of 7 versus a with 6 — 0.4, c\ — 2, c2 — 2.1 (Fig. 5). Clearly two curves of Hopf bifurcation persist, but they no longer intersect. The picture is similar for other values of the parameters, with the distance between the curves being a function of the difference between c\ and c2. A measure of this distance (at 7 — 0) is ajji — ajj2 where the ajjj are given by (33).
(a)
0.78	
0.62	
(b) 1.5
-1.5
0.14
-0.1
Figure 5 c2 =2.1, fe^igram
1.5
Yi
0.1
, (a) Hopf bifurcation sets in 7, a plane for non-identical oscillators (30). Here b — 0.4, c\ — 2, but the picture would be qualitatively similar for any b < min(l, \cj\), c\ ^ c^- (b) Bifurcation of Xi vs 7 for the same parameter values ari&)a = 0.7.
1.5
Figure 6. Numerical simulations of (30) with a — 0.7, b — 0.4, c\ — 2 and — 2.1 for three initial conditions resulting in three different attractors. (a) Attractors are all limit cycles for 7 = 0.05 (b) Two attractors have undergone period doubling bifurcations at 7 = 0.076.
As discussed above, changing the value of c's in the model is the "least harmful" way of breaking the symmetry of the equations. Since this resulted in a loss of Hopf interaction points we might expect the same to occur for other changes to the parameters. Indeed, this is what we observe when we change either the 6's or the a's.
Appealing to the continuous dependence on parameters of solutions of well behaved ODE's, we might expect that behaviour in the symmetric system would persist in the non-symmetric system if the parameters of the two neurons are close. In fact, for the parameter values near to the "closest approach" of the two Hopf curves of Fig. 5 we observe the coexistence of three limit cycles in the systems. For slightly large values of the coupling parameter 7 we observe the coexistence of a large amplitude limit cycles with two smaller amplitude high period attractors (Fig. 6). For a larger difference between the cj's we observe bistability between two limit cycles or a limit cycle and a more complicated attractor or two complicated attractors. This behaviour exists for small values of the coupling even if the the difference between the two parameters is up to 20% of the size of the parameter. We observe similar behaviour when the a's or the 6's are chosen to differ and when multiple parameters differ.
5. Discussion
There have been many studies of coupled relaxation oscillators, including the Fitzhugh-Nagumo [6] and van der Pol [8-15] models (see also references therein). Much of this work concerns the existence
10
and stability of various types (in phase, out of phase, phase locked) of periodic motions, often in the case identical oscillators. By contrast, our work focussed primarily of the origin of multistability in a model with identical neurons and whether this multistability persists when the neurons are no longer identical.
We have shown that the coexistence of multiple stable limit cycles or stable limit cycles and other more complicated attractors result from a symmetric Hopf bifurcation when two identical neurons are coupled together. Similar behaviour has been observed in nonlinearly coupled van der Pol oscillators [15] and Fitzhugh-Nagumo equations with electrical coupling in both the Xj and yj [6]. In fact such behaviour can be expected to occur for any model which can undergo a Hopf bifurcation [16].
When the neural models are identical, but the coupling is no longer symmetric we have shown that the symmetric Hopf bifurcation becomes a one-to-one resonant Hopf-Hopf interaction, and multistability still occurs. This result also only depends on having two identical neurons which can undergo a Hopf bifurcation, and thus should also hold for other neural models.
Finally we showed that when the slightest change in one of the parameters in the neural models is introduced the Hopf interaction point is totally lost. However, if the difference of the parameters is no too great the multistability remains. We conjecture that this result will also hold for arbitrary neural models with electrical coupling.
We have not emphasized the fact, but it should be clear that the multistability we have observed occurs close to the interaction points. Generally, this means that the coupling between the neurons must be small to observe such behaviour. This is in agreement with the work in [?].
6. Acknowledgements
The numerical simulations Figs. 2, 4, 6 were performed with XPPAUT [17]. The bifurcation diagrams for Figs. 1(b), 3, 4 were created using AUTO [18] through the XPPAUT interface. SAC would like to thank the Centre de recherches mathematiques (Universite de Montreal) for hospitality during the writing of this paper.
REFERENCES
1. R. Fitzhugh, Thresholds and plateaus in the Hodgkin-Huxley nerve equations, Journal of General Physiology 43 (1960) 867.
2. J. Nagumo, S. Arimoto, S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proc. IRE 50 (1962) 2061-2070.
3. B. van der Pol, On relaxation-oscillations, Phil. Mag. 2 (1926) 978-992, 7th Series.
4. A. Hodgkin, A. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, Journal of Physiology 117 (1952) 500-544.
5. F. Skinner, L. Zhang, J. Velazquez, P. Carlen, Bursting in inhibitory interneuronal networks: A role for gap-junctional coupling, Journal of Neurophysiology 81 (3) (1999) 1274-83.
6. M. Kawato, M. Sokabe, R. Suzuki, Synergism and antagonism of neurons caused by an electrical synapse, Biological Cybernetics 34 (1979) 81-89.
7. M. Golubitsky, I. Stewart, D. Schaeffer, Singularities and Groups in Bifurcation Theory, Vol. 2, Springer Verlag, New York, 1988.
8. V. Arnol'd, Geometrical Methods in the Theory of Ordinary Differential Equations, Springer Verlag, New York, 1983.
9. J. K. Aggarwal, C. G. Richie, On coupled van der Pol oscillators, IEEE Trans. Circuit Theory CT-13 (1966) 465-466.
10. J. Belair, P. Holmes, On linearly coupled relaxation oscillators, Quarterly of Applied Mathematics February (1984) 193-219.
11. T. Chakraborty, R. H. Rand, The transition from phase locking to drift in a system of two weakly coupled van der Pol oscillators, Internat. J. Non-Linear Mech. 23 (5-6) (1988) 369-376.
12. J. Grasman, M. J. W. Jansen, Mutually synchronized relation oscillators as prototypes of oscillating systems in biology, J. Math. Biol. 7 (2) (1979) 171-197.
13. T. Kawahara, Coupled van der Pol oscillators—a model of excitatory and inhibitory neural interactions, Biol. Cybernet. 39 (1) (1980/81) 37-43.
14. A. K. Kozlov, M. M. Sushchik, Y. I. Molkov, A. S. Kuznetsov, Bistable phase synchronization and
11
chaos in a system of coupled van der Pol-Duffing oscillators, Internat. J. Bifur. Chaos Appl. Sei. Engrg. 9 (12) (1999) 2271-2277.
15. E. Palm, M. Tveitereid, On coupled van der Pol equations, Quart. J. Mech. Appl. Math. 33 (3) (1980) 267-276.
16. I. Pastor-Diaz, A. Lopez-Fraguas, Dynamics of two coupled van der Pol oscillators, Physical Review E 52 (2) (1995) 1480-1489.
17. F. Hoppensteadt, E. Izhikevich, Weakly connected neural networks, Springer-Verlag, New York, 1997.
18. B. Ermentrout, XPPAUT3.0 - the differential equations tool, Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, 1997.
19. E. Doedel, Auto: A program for the automatic bifurcation analysis of autonomous systems, in: Proc. 10th Manitoba Conf. on Num. Math, and Comp., Univ. of Manitoba, Winnipeg, Canada, 1981, pp. 265-284.