Table Of ContentHomoclinic Orbits of the FitzHugh-Nagumo
Equation: The Singular-Limit
2
1
0
2 John Guckenheimer and Christian Kuehn
n
a
J Abstract
7
2 TheFitzHugh-Nagumoequationhasbeeninvestigatedwithawidear-
ray of different methods in the last three decades. Recently a version of
] theequationswithanappliedcurrentwasanalyzedbyChampneys,Kirk,
S
Knobloch,OldemanandSneyd[5]usingnumericalcontinuationmethods.
D
Theyobtainedacomplicatedbifurcationdiagraminparameterspacefea-
. turingaC-shapedcurveofhomoclinic bifurcationsandaU-shapedcurve
h
ofHopfbifurcations. Weusetechniquesfrom multipletime-scale dynam-
t
a icstounderstandthestructuresofthisbifurcationdiagrambasedongeo-
m
metric singular perturbation analysis of the FitzHugh-Nagumo equation.
[ Numerical and analytical techniques show that if the ratio of the time-
scales in theFitzHugh-Nagumoequation tendstozero, thenoursingular
1
limit analysis correctly represents the observed CU-structure. Geometric
v
insight from theanalysis can evenbeused to computebifurcation curves
1
which areinaccessible viacontinuationmethods. Theresultsof ouranal-
0
ysis are summarized in a singular bifurcation diagram.
9
5
.
1 1 Introduction
0
2
1.1 Fast-Slow Systems
1
:
v Fast-slow systems of ordinary differential equations (ODEs) have the general
i form:
X
r dx
a ǫx˙ = ǫ =f(x,y,ǫ) (1)
dτ
dy
y˙ = =g(x,y,ǫ)
dτ
where x Rm, y Rn and 0 ǫ 1 represents the ratio of time scales. The
∈ ∈ ≤ ≪
functions f and g are assumed to be sufficiently smooth. In the singular limit
ǫ 0 the vector field (1) becomes a differential-algebraic equation. The alge-
br→aic constraint f = 0 defines the critical manifold C = (x,y) Rm Rn :
0
{ ∈ ×
f(x,y,0) = 0 . Where D f(p) is nonsingular, the implicit function theorem
x
}
implies that there exists a map h(x) = y parametrizing C as a graph. This
0
yields the implicitly defined vector field y˙ = g(h(y),y,0) on C called the slow
0
1
flow.
We canchange(1) to the fasttime scale t=τ/ǫ andlet ǫ 0 to obtainthe
→
second possible singular limit system
dx
x = =f(x,y,0) (2)
′
dt
dy
y = =0
′
dt
Wecallthevectorfield(2)parametrizedbytheslowvariablesy thefastsubsys-
tem or the layer equations. The central idea of singular perturbation analysis
is to use information about the fast subsystem and the slow flow to under-
stand the full system (1). One of the main tools is Fenichel’s Theorem (see
[16, 17, 18, 19]). It states that for every ǫ sufficiently small and C normally
0
hyperbolic there exists a family of invariant manifolds C for the flow (1). The
ǫ
manifoldsareatadistanceO(ǫ)fromC andthe flowsonthemconvergetothe
0
slow flow on C as ǫ 0. Points p C where D f(p) is singular are referred
0 0 x
→ ∈
to as fold points1.
Beyond Fenichel’s Theorem many other techniques have been developed.
More detailed introductions and results can be found in [12, 34, 24] from a
geometric viewpoint. Asymptotic methods are developed in [42, 23] whereas
ideas from nonstandard analysis are introduced in [8]. While the theory is well
developed for two-dimensional fast-slow systems, higher-dimensional fast-slow
systemsareanactiveareaofcurrentresearch. Inthefollowingweshallfocuson
theFitzHugh-Nagumoequationviewedasathree-dimensionalfast-slowsystem.
1.2 The FitzHugh-Nagumo Equation
TheFitzHugh-NagumoequationisasimplificationoftheHodgin-Huxleymodel
for an electric potential of a nerve axon [30]. The first version was developed
by FitzHugh [20] and is a two-dimensional system of ODEs:
u3
ǫu˙ = v +u+p (3)
− 3
1
v˙ = (v+γu a)
−s −
A detailed summary of the bifurcations of (3) can be found in [44]. Nagumo et
al. [43]studiedarelatedequationthatadds adiffusiontermforthe conduction
process of action potentials along nerves:
u =δu +f (u) w+p
τ xx a − (4)
w =ǫ(u γw)
τ
(cid:26) −
1TheprojectionofC0 ontothexcoordinatesmayhavemoredegeneratesingularitiesthan
foldsingularitiesatsomeofthesepoints.
2
where f (u) = u(u a)(1 u) and p,γ,δ and a are parameters. A good in-
a
− −
troduction to the derivation and problems associated with (4) can be found in
[28]. Suppose we assume a traveling wave solution to (4) and set u(x,τ) =
u(x+sτ)=u(t) and w(x,τ) =w(x+sτ) =w(t), where s represents the wave
speed. By the chain rule we get u = su, u = u and w = sw. Set v = u
τ ′ xx ′′ τ ′ ′
and substitute into (4) to obtain the system:
u = v
′
1
v′ = (sv fa(u)+w p) (5)
δ − −
ǫ
w = (u γw)
′
s −
System (5) is the FitzHugh-Nagumo equation studied in this paper. Observe
that a homoclinic orbit of (5) corresponds to a traveling pulse solution of (4).
These solutions are of special importance in neuroscience [28] and have been
analyzedusing severaldifferent methods. For example, it has been provedthat
(5) admits homoclinic orbits [29, 4] for small wave speeds (“slow waves”) and
largewavespeeds (“fast waves”). Fastwavesare stable [33] andslow wavesare
unstable [21]. It has been shown that double-pulse homoclinic orbits [15] are
possible. If (5) has two equilibrium points and heteroclinic connections exist,
bifurcation from a twisted double heteroclinic connection implies the existence
ofmulti-pulsetravelingfrontandbackwaves[6]. Theseresultsarebasedonthe
assumptionofcertainparameterrangesforwhichwerefertotheoriginalpapers.
Geometric singular perturbation theory has been used successfully to analyze
(5). In [32] the fast pulse is constructed using the exchange lemma [35, 31, 3].
The exchange lemma has also been used to prove the existence of a codimen-
sion two connection between fast and slow waves in (s,ǫ,a)-parameter space
[38]. An extension of Fenichel’s theorem and Melnikov’s method can be em-
ployed to prove the existence of heteroclinic connections for parameter regimes
of (5) with two fixed points [45]. The general theory of relaxation oscillations
in fast-slow systems applies to (5) (see e.g. [42, 26]) as does - at least partially
- the theory of canards (see e.g. [46, 9, 11, 39]).
The equations (5) have been analyzed numerically by Champneys, Kirk,
Knobloch, Oldeman and Sneyd [5] using the numerical bifurcation software
AUTO [13, 14]. They considered the following parameter values:
1
γ =1, a= , δ =5
10
We shall fix those values to allow comparison of our results with theirs. Hence
we also write f (u) = f(u). Changing from the fast time t to the slow time
1/10
τ and relabeling variables x =u, x =v and y =w we get:
1 2
ǫx˙ = x
1 2
1 1 1
ǫx˙ = (sx x (x 1)( x )+y p)= (sx f(x )+y p)(6)
2 2 1 1 1 2 1
5 − − 10 − − 5 − −
1
y˙ = (x y)
1
s −
3
Fromnowonwereferto(6)as“the”FitzHugh-Nagumoequation. Investigating
bifurcations in the (p,s) parameter space one finds C-shaped curves of homo-
clinic orbits and a U-shaped curve of Hopf bifurcations; see Figure 1. Only
part of the bifurcation diagram is shown in Figure 1. There is another curve of
homoclinic bifurcations onthe rightside of the U-shaped Hopf curve. Since (6)
has the symmetry
11 11 11 33
x x , x x , y y, p 1 p (7)
1 1 2 2
→ 15 − → 15 − →− → 15 − 225 −
(cid:18) (cid:19)
we shall examine only the left side of the U-curve. The homoclinic C-curve is
difficult to compute numerically by continuation methods using AUTO [13, 14]
orMatCont[22]. Thecomputationsseeminfeasibleforsmallvaluesofǫ 10 3.
−
≤
Furthermore multi-pulse homoclinic orbits can exist very close to single pulse
onesanddistinguishingbetweenthemmustnecessarilyencounterproblemswith
numerical precision [5]. The Hopf curve and the bifurcations of limit cycles
showninFigure1havebeencomputedusingMatCont. Thecurveofhomoclinic
bifurcationshasbeencomputedbyanewmethodtobedescribedinSection3.2.
2.2
2
Hopf ε=0.01
1.8
1.6
SNLC
1.4
s
1.2
1
Hom
0.8
0.6 GH
0.4
0 0.05 0.1 p 0.15 0.2 0.25
Figure 1: Bifurcation diagram of (6). Hopf bifurcations are shown in green,
saddle-node of limit cycles (SNLC) are shown in blue and GH indicates a gen-
eralized Hopf (or Bautin) bifurcation. The arrows indicate the side on which
periodic orbits are generated at the Hopf bifurcation. The red curve shows
(possible) homoclinic orbits; in fact, homoclinic orbits only exist to the left of
the twoblackdots(see Section3.2). Only partofthe parameterspaceis shown
because ofthe symmetry (7). The homoclinic curve hasbeen thickenedto indi-
cate that multipulse homoclinic orbits exist very close to single pulse ones (see
[15]).
SincethebifurcationstructureshowninFigure1wasalsoobservedforother
4
excitable systems,Champneys et al.[5] introducedthe term CU-system. Bifur-
cationanalysisfromtheviewpointofgeometricsingularperturbationtheoryhas
beencarriedoutforexampleswithonefastandtwoslowvariables[27,2,25,41].
Since the FitzHugh-Nagumo equation has one slow and two fast variables, the
situation is quite different and new techniques have to be developed. Our main
goal is to show that many features of the complicated 2-parameter bifurcation
diagramshowninFigure1canbederivedwithacombinationoftechniquesfrom
singularperturbationtheory,bifurcationtheoryandrobustnumericalmethods.
We accurately locate where the system has canards and determine the orbit
structure of the homoclinic and periodic orbits associated to the C-shaped and
U-shaped bifurcation curves, without computing the canards themselves. We
demonstrate that the basic CU-structure of the system can be computed with
elementarymethodsthatdonotusecontinuationmethodsbasedoncollocation.
The analysis of the slow and fast subsystems yields a “singularbifurcation dia-
gram” to which the basic CU structure in Figure 1 converges as ǫ 0.
→
Remark: We have also investigated the termination mechanism of the C-
shaped homoclinic curve described in [5]. Champneys et al. observed that the
homocliniccurvedoesnotreachtheU-shapedHopfcurvebutturnsaroundand
foldsbackclosetoitself. Wecomputeaccurateapproximationsofthehomoclinic
orbits for smaller values ǫ than seems possible with AUTO in this region. One
aspectofouranalysisisanewalgorithmforcomputinginvariantslowmanifolds
of saddle type in the full system. This work will be described elsewhere.
2 The Singular Limit
The first step in our analysis is to investigate the slow and fast subsystems
separately. Let ǫ 0 in (6); this yields two algebraic constraints that define
→
the critical manifold:
1
C = (x ,x ,y) R3 :x =0 y =x (x 1)( x )+p=c(x )
0 1 2 2 1 1 1 1
∈ − 10 −
(cid:26) (cid:27)
Therefore C is a cubic curve in the coordinate plane x =0. The parameter p
0 2
moves the cubic up and down inside this plane. The critical points of the cubic
are solutions of c(x )=0 and are given by:
′ 1
1
x = 11 √91 or numerically: x 0.6846, x 0.0487
1, 1,+ 1,
± 30 ± ≈ − ≈
(cid:16) (cid:17)
Thepointsx arefoldpointswith c (x ) =0sinceC isacubicpolynomial
1, ′′ 1, 0
± | ± |6
with distinct critical points. The fold points divide C into three segments
0
C = x <x C , C = x x x C , C = x <x C
l 1 1, 0 m 1, 1 1,+ 0 r 1,+ 1 0
{ −}∩ { − ≤ ≤ }∩ { }∩
We denote the associated slow manifolds by C , C and C . There are
l,ǫ m,ǫ r,ǫ
two possibilities to obtain the slow flow. One way is to solve c(x ) = y for
1
5
x and substitute the result into the equation y˙ = 1(x y). Alternatively
1 s 1 −
differentiating y = c(x ) implicitly with respect to τ yields y˙ = x˙ c(x ) and
1 1 ′ 1
therefore
1 1
(x y)=x˙ c(x ) x˙ = (x c(x )) (8)
1 1 ′ 1 1 1 1
s − ⇒ sc(x ) −
′ 1
One can view this as a projection of the slow flow, which is constrained to the
critical manifold in R3, onto the x -axis. Observe that the slow flow is singu-
1
lar at the fold points. Direct computation shows that the fixed point problem
x = c(x ) has only a single real solution. This implies that the critical mani-
1 1
fold intersects the diagonaly =x only in a single point x which is the unique
1 ∗1
equilibrium of the slow flow (8). Observe that q =(x ,0,x ) is also the unique
∗1 ∗1
equilibrium of the full system (6) and depends on p. Increasing p moves the
equilibrium from left to right on the critical manifold. The easiest practical
way to determine the direction of the slow flow on C is to look at the sign of
0
(x y). The situation is illustrated in Figure 2.
1
−
Slow Flow in the plane x=0
2
0.4
y=x
1
0.3
0.2
eq. pt. q
0.1 x
1,+
y 0
−0.1
x
1,−
−0.2
C
−0.3 0
−0.4
−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2
x
1
Figure 2: Sketch of the slow flow on the critical manifold C
0
2.1 The Slow Flow
Weareinterestedinthebifurcationsoftheslowflowdependingontheparameter
p. The bifurcations occur when x passes through the fold points. The values
∗1
ofp cansimply be foundby solvingthe equationsc(x )=0 andc(x ) x =0
′ 1 1 1
−
simultaneously. The result is:
p 0.0511 and p 0.5584
+
− ≈ ≈
wherethesubscriptsindicatethefoldpointatwhicheachequilibriumislocated.
6
The singular time-rescaling τ¯ = sc(x )/τ of the slow flow yields the desin-
′ 1
gularized slow flow
dx x
1 1
=x c(x )=x + (x 1)(10x 1) p (9)
1 1 1 1 1
dτ¯ − 10 − − −
Time is reversed by this rescaling on C and C since s > 0 and c(x ) is
l r ′ 1
negativeonthese branches. The desingularizedslowflow (9)is smoothandhas
no bifurcations as p is varied.
2.2 The Fast Subsystem
Thekeycomponentofthefast-slowanalysisfortheFitzHugh-Nagumoequation
is the two-dimensional fast subsystem
x = x
′1 2
1 1
x = (sx x (x 1)( x )+y p) (10)
′2 5 2− 1 1− 10 − 1 −
where p 0, s 0 are parameters and y is fixed. Since y and p have the same
≥ ≥
effect as bifurcation parameters we set p y = p¯. We consider several fixed
−
y-values and the effect of varying p (cf. Section 3.2) in each case. There are
either one, two or three equilibrium points for (10). Equilibrium points satisfy
x =0 and lie on the critical manifold, i.e. we have to solve
2
1
0=x (x 1)( x )+p¯ (11)
1 1 1
− 10 −
forx . Wefindthattherearethreeequilibriaforapproximatelyp¯ = 0.1262<
1 l
−
p¯ < 0.0024 = p¯ , two equilibria on the boundary of this p interval and one
r
equilibrium otherwise. The Jacobian of (10) at an equilibrium is
0 1
A(x )=
1 1 1 22x +30x2 s
(cid:18) 50 − 1 1 5 (cid:19)
(cid:0) (cid:1)
Direct calculation yields that for p [p¯,p¯ ] the single equilibrium is a saddle.
l r
6∈
In the case of three equilibria, we have a source that lies between two saddles.
Note that this also describes the stability of the three branches of the critical
manifoldC ,C andC . Fors>0thematrixAissingularofrank1ifandonly
l m r
if 30x2 22x +1=0 which occurs for the fold points x . Hence the equilib-
1− 1 1,±
ria of the fast subsystem undergo a fold (or saddle-node) bifurcation once they
approach the fold points of the critical manifold. This happens for parameter
valuesp¯ andp¯ . Notethatbysymmetrywecanreducetostudyingasinglefold
l r
point. In the limit s = 0 (corresponding to the case of a “standing wave”) the
saddle-node bifurcation point becomes more degenerate with A(x ) nilpotent.
1
7
Our next goal is to investigate global bifurcations of (10); we start with
homoclinic orbits. For s=0 itis easyto see that (10) is a Hamiltoniansystem:
∂H
x = =x
′1 ∂x 2
2
∂H 1 1
x = = ( x (x 1)( x ) p¯) (12)
′2 −∂x 5 − 1 1− 10 − 1 −
1
with Hamiltonian function
1 (x )2 11(x )3 (x )4 x p¯
H(x ,x )= x2 1 + 1 1 + 1 (13)
1 2 2 2− 100 150 − 20 5
We will use this Hamiltonian formulation later on to describe the geometry of
homoclinic orbits for slow wave speeds. Assume that p¯ is chosen so that (12)
has a homoclinic orbit x (t). We are interested in perturbations with s > 0
0
and note that in this case the divergence of (10) is s. Hence the vector field
is area expanding everywhere. The homoclinic orbit breaks for s > 0 and no
periodic orbits are created. Note that this scenario does not apply to the full
three-dimensional system as the equilibrium q has a pair of complex conjugate
eigenvalues so that a Shilnikov scenario can occur. This illustrates that the
singular limit can be used to help locate homoclinic orbits of the full system,
but that some characteristics of these orbits change in the singular limit.
We are interested next in finding curves in (p¯,s)-parameter space that rep-
resent heteroclinic connections of the fast subsystem. The main motivation is
thedecompositionoftrajectoriesinthe fullsystemintoslowandfastsegments.
Concatenating fast heteroclinic segments and slow flow segments can yield ho-
moclinicorbitsofthefullsystem[28,4,32,38]. Wedescribeanumericalstrategy
to detect heteroclinic connections in the fast subsystem and continue them in
parameter space. Suppose that p¯ (p¯,p¯ ) so that (10) has three hyperbolic
l r
∈
equilibrium points x, x and x . We denote by Wu(x ) the unstable and by
l m r l
Ws(x ) the stable manifold of x . The same notation is also used for x and
l l r
tangentspacesto Ws(.)andWu(.)aredenotedby Ts(.)andTu(.). Recallthat
x isasourceandshallnotbeofinteresttousfornow. Definethecrosssection
m
Σ by
x +x
Σ= (x ,x ) R2 :x = l r .
1 2 1
{ ∈ 2 }
WeuseforwardintegrationofinitialconditionsinTu(x )andbackwardintegra-
l
tionofinitialconditionsinTs(x )toobtaintrajectoriesγ+ andγ respectively.
r −
We calculate their intersection with Σ and define
γl(p¯,s):=γ+ Σ, γr(p¯,s):=γ− Σ
∩ ∩
We compute the functions γ and γ for different parameter values of (p¯,s)
l r
numerically. Heteroclinic connections occur at zeros of the function
h(p¯,s):=γ (p¯,s) γ (p¯,s)
l r
−
8
Oncewefindaparameterpair(p¯ ,s )suchthath(p¯ ,s )=0,theseparameters
0 0 0 0
can be continued along a curve of heteroclinic connections in (p¯,s) parameter
space by solving the root-finding problem h(p¯ +δ ,s +δ ) = 0 for either δ
0 1 0 2 1
or δ fixed and small. We use this method later for different fixed values of y
2
to compute heteroclinic connections in the fast subsystem in (p,s) parameter
space. The results of these computations are illustrated in Figure 3. There are
two distinct branches in Figure 3. The branches are asymptotic to p¯ and p¯
l r
and approximately form a “V”. From Figure 3 we conjecture that there exists
a double heteroclinic orbit for p¯ 0.0622.
≈−
1.5
1
s
0.5
0
−0.15 −0.1 −0.05 0 0.05
p¯
Figure 3: Heteroclinic connections for equation (10) in parameter space.
Remarks: If we fix p=0 our initial change of variable becomes y =p¯and
−
our results for heteroclinic connections are for the FitzHugh-Nagumo equation
without an applied current. In this situation it has been shown that the hete-
roclinic connections of the fast subsystem can be used to prove the existence of
homoclinic orbits to the unique saddle equilibrium (0,0,0)(cf. [32]). Note that
the existence of the heteroclinics in the fast subsystem was proved in a special
caseanalytically[1]but Figure3 is - to the best ofourknowledge- the firstex-
plicitcomputationofwherefastsubsystemheteroclinicsarelocated. Thepaper
[36] develops a method for finding heteroclinic connections by the same basic
approachwe used, i.e. defining acodimensionone hyperplane H that separates
equilibrium points.
Figure3suggeststhatthereexistsadoubleheteroclinicconnectionfors=0.
Observe that the Hamiltonian in our case is H(x ,x ) = (x2)2 +V(x ) where
1 2 2 1
the function V(x ) is:
1
px (x )2 11(x )3 (x )4
1 1 1 1
V(x )= +
1
5 − 100 150 − 20
The solution curves of (12) are given by x = 2(const. V(x )). The
2 1
± −
structureofthe solutioncurvesentailssymmetryunderreflectionaboutthe x -
p 1
axis. Suppose p¯ [p¯,p¯ ] and recall that we denoted the two saddle points of
l r
∈
(10) by x andx andthat their locationdepends onp¯. Therefore,we conclude
l r
9
that the two saddles x and x must have a heteroclinic connection if they lie
l r
on the same energy level, i.e. they satisfy V(x ) V(x ) = 0. This equation
l r
−
can be solved numerically to very high accuracy.
Proposition 1. The fast subsystem of the FitzHugh-Nagumo equation for s=
0 has a double heteroclinic connection for p¯ = p¯ 0.0619259. Given a
∗
≈ −
particular value y = y there exists a double heteroclinic connection for p =
0
p¯ +y in the fast subsystem lying in the plane y =y .
∗ 0 0
2.3 Two Slow Variables, One Fast Variable
From continuation of periodic orbits in the full system - to be described in
Section 3.1 - we observe that near the U-shaped curve of Hopf bifurcations the
x -coordinate is a faster variable than x . In particular, the small periodic
2 1
orbits generated in the Hopf bifurcation lie almost in the plane x = 0. Hence
2
to analyze this region we set x¯ = x /ǫ to transform the FitzHugh-Nagumo
2 2
equation (6) into a system with 2 slow and 1 fast variable:
x˙ = x¯
1 2
1 1
ǫ2x¯˙ = (sǫx¯ x (x 1)( x )+y p) (14)
2 2 1 1 1
5 − − 10 − −
1
y˙ = (x y)
1
s −
Note that (14) corresponds to the FitzHugh-Nagumo equation in the form (cf.
(4)):
u =5ǫ2u +f(u) w+p
τ xx − (15)
w =ǫ(u w)
τ
(cid:26) −
Therefore the transformation x¯ = x /ǫ can be viewed as a rescaling of the
2 2
diffusion strength by ǫ2. We introduce a new independent small parameter
δ¯=ǫ2 and then let δ¯=ǫ2 0. This assumes that O(ǫ) terms do not vanish in
→
this limit, yielding the diffusion free system. Then the slow manifold S of (14)
0
is:
1
S = (x ,x¯ ,y) R3 :x¯ = (f(x ) y+p) (16)
0 1 2 2 1
∈ sǫ −
(cid:26) (cid:27)
Proposition 2. Following time rescaling by s, the slow flow of system (14) on
S in the variables (x ,y) is given by
0 1
ǫx˙ = f(x ) y+p
1 1
−
y˙ = x y (17)
1
−
In the variables (x ,x¯ ) the vector field (17) becomes
1 2
x˙ = x¯
1 2
1 x¯
ǫx¯˙ = (x f(x ) p)+ 2 (f (x ) ǫ) (18)
2 −s2 1− 1 − s ′ 1 −
10
