Homoclinic 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