Homoclinic Orbits of the FitzHugh-Nagumo Equation: The Singular-Limit

Homoclinic Orbits of the FitzHugh-Nagumo Equation: The Singular-Limit

John Guckenheimer and Christian Kuehn
Abstract

The FitzHugh-Nagumo equation has been investigated with a wide array of different methods in the last three decades. Recently a version of the equations with an applied current was analyzed by Champneys, Kirk, Knobloch, Oldeman and Sneyd [5] using numerical continuation methods. They obtained a complicated bifurcation diagram in parameter space featuring a C-shaped curve of homoclinic bifurcations and a U-shaped curve of Hopf bifurcations. We use techniques from multiple time-scale dynamics to understand the structures of this bifurcation diagram based on geometric singular perturbation analysis of the FitzHugh-Nagumo equation. Numerical and analytical techniques show that if the ratio of the time-scales in the FitzHugh-Nagumo equation tends to zero, then our singular limit analysis correctly represents the observed CU-structure. Geometric insight from the analysis can even be used to compute bifurcation curves which are inaccessible via continuation methods. The results of our analysis are summarized in a singular bifurcation diagram.

1 Introduction

1.1 Fast-Slow Systems

Fast-slow systems of ordinary differential equations (ODEs) have the general form:

 ϵ˙x = ϵdxdτ=f(x,y,ϵ) (1) ˙y = dydτ=g(x,y,ϵ)

where , and represents the ratio of time scales. The functions and are assumed to be sufficiently smooth. In the singular limit the vector field (1) becomes a differential-algebraic equation. The algebraic constraint defines the critical manifold . Where is nonsingular, the implicit function theorem implies that there exists a map parametrizing as a graph. This yields the implicitly defined vector field on called the slow flow.

We can change (1) to the fast time scale and let to obtain the second possible singular limit system

 x′ = dxdt=f(x,y,0) (2) y′ = dydt=0

We call the vector field (2) parametrized by the slow variables the fast subsystem or the layer equations. The central idea of singular perturbation analysis is to use information about the fast subsystem and the slow flow to understand 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 normally hyperbolic there exists a family of invariant manifolds for the flow (1). The manifolds are at a distance from and the flows on them converge to the slow flow on as . Points where is singular are referred to as fold points111The projection of onto the coordinates may have more degenerate singularities than fold singularities at some of these points..

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 systems are an active area of current research. In the following we shall focus on the FitzHugh-Nagumo equation viewed as a three-dimensional fast-slow system.

1.2 The FitzHugh-Nagumo Equation

The FitzHugh-Nagumo equation is a simplification of the Hodgin-Huxley model 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:

 ϵ˙u = v−u33+u+p (3) ˙v = −1s(v+γu−a)

A detailed summary of the bifurcations of (3) can be found in [44]. Nagumo et al. [43] studied a related equation that adds a diffusion term for the conduction process of action potentials along nerves:

 {uτ=δuxx+fa(u)−w+pwτ=ϵ(u−γw) (4)

where and and are parameters. A good introduction to the derivation and problems associated with (4) can be found in [28]. Suppose we assume a traveling wave solution to (4) and set and , where represents the wave speed. By the chain rule we get , and . Set and substitute into (4) to obtain the system:

 u′ = v v′ = 1δ(sv−fa(u)+w−p) (5) w′ = ϵs(u−γw)

System (1.2) is the FitzHugh-Nagumo equation studied in this paper. Observe that a homoclinic orbit of (1.2) corresponds to a traveling pulse solution of (4). These solutions are of special importance in neuroscience [28] and have been analyzed using several different methods. For example, it has been proved that (1.2) admits homoclinic orbits [29, 4] for small wave speeds (“slow waves”) and large wave speeds (“fast waves”). Fast waves are stable [33] and slow waves are unstable [21]. It has been shown that double-pulse homoclinic orbits [15] are possible. If (1.2) has two equilibrium points and heteroclinic connections exist, bifurcation from a twisted double heteroclinic connection implies the existence of multi-pulse traveling front and back waves [6]. These results are based on the assumption of certain parameter ranges for which we refer to the original papers. Geometric singular perturbation theory has been used successfully to analyze (1.2). 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 codimension two connection between fast and slow waves in ()-parameter space [38]. An extension of Fenichel’s theorem and Melnikov’s method can be employed to prove the existence of heteroclinic connections for parameter regimes of (1.2) with two fixed points [45]. The general theory of relaxation oscillations in fast-slow systems applies to (1.2) (see e.g. [42, 26]) as does - at least partially - the theory of canards (see e.g. [46, 9, 11, 39]).

The equations (1.2) 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,a=110,δ=5

We shall fix those values to allow comparison of our results with theirs. Hence we also write . Changing from the fast time to the slow time and relabeling variables , and we get:

 ϵ˙x1 = x2 ϵ˙x2 = 15(sx2−x1(x1−1)(110−x1)+y−p)=15(sx2−f(x1)+y−p) (6) ˙y = 1s(x1−y)

From now on we refer to (1.2) as “the” FitzHugh-Nagumo equation. Investigating bifurcations in the () parameter space one finds C-shaped curves of homoclinic 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 on the right side of the U-shaped Hopf curve. Since (1.2) has the symmetry

 x1→1115−x1,x2→1115−x2,y→−y,p→1115(1−33225)−p (7)

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] or MatCont [22]. The computations seem infeasible for small values of . Furthermore multi-pulse homoclinic orbits can exist very close to single pulse ones and distinguishing between them must necessarily encounter problems with numerical precision [5]. The Hopf curve and the bifurcations of limit cycles shown in Figure 1 have been computed using MatCont. The curve of homoclinic bifurcations has been computed by a new method to be described in Section 3.2.

Since the bifurcation structure shown in Figure 1 was also observed for other excitable systems, Champneys et al. [5] introduced the term CU-system. Bifurcation analysis from the viewpoint of geometric singular perturbation theory has been carried out for examples with one fast and two slow variables [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 diagram shown in Figure 1 can be derived with a combination of techniques from singular perturbation theory, bifurcation theory and robust numerical methods. 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 elementary methods that do not use continuation methods based on collocation. The analysis of the slow and fast subsystems yields a “singular bifurcation diagram” to which the basic CU structure in Figure 1 converges as .

Remark: We have also investigated the termination mechanism of the C-shaped homoclinic curve described in [5]. Champneys et al. observed that the homoclinic curve does not reach the U-shaped Hopf curve but turns around and folds back close to itself. We compute accurate approximations of the homoclinic orbits for smaller values than seems possible with AUTO in this region. One aspect of our analysis is a new algorithm for computing invariant slow manifolds 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 in (1.2); this yields two algebraic constraints that define the critical manifold:

 C0={(x1,x2,y)∈R3:x2=0y=x1(x1−1)(110−x1)+p=c(x1)}

Therefore is a cubic curve in the coordinate plane . The parameter moves the cubic up and down inside this plane. The critical points of the cubic are solutions of and are given by:

 x1,±=130(11±√91)or numerically:x1,+≈0.6846,x1,−≈0.0487

The points are fold points with since is a cubic polynomial with distinct critical points. The fold points divide into three segments

 Cl={x1

We denote the associated slow manifolds by , and . There are two possibilities to obtain the slow flow. One way is to solve for and substitute the result into the equation . Alternatively differentiating implicitly with respect to yields and therefore

 1s(x1−y)=˙x1c′(x1)⇒˙x1=1sc′(x1)(x1−c(x1)) (8)

One can view this as a projection of the slow flow, which is constrained to the critical manifold in , onto the -axis. Observe that the slow flow is singular at the fold points. Direct computation shows that the fixed point problem has only a single real solution. This implies that the critical manifold intersects the diagonal only in a single point which is the unique equilibrium of the slow flow (8). Observe that is also the unique equilibrium of the full system (1.2) and depends on . Increasing moves the equilibrium from left to right on the critical manifold. The easiest practical way to determine the direction of the slow flow on is to look at the sign of . The situation is illustrated in Figure 2.

2.1 The Slow Flow

We are interested in the bifurcations of the slow flow depending on the parameter . The bifurcations occur when passes through the fold points. The values of can simply be found by solving the equations and simultaneously. The result is:

 p−≈0.0511andp+≈0.5584

where the subscripts indicate the fold point at which each equilibrium is located.

The singular time-rescaling of the slow flow yields the desingularized slow flow

 dx1d¯τ=x1−c(x1)=x1+x110(x1−1)(10x1−1)−p (9)

Time is reversed by this rescaling on and since and is negative on these branches. The desingularized slow flow (9) is smooth and has no bifurcations as is varied.

2.2 The Fast Subsystem

The key component of the fast-slow analysis for the FitzHugh-Nagumo equation is the two-dimensional fast subsystem

 x′1 = x2 x′2 = 15(sx2−x1(x1−1)(110−x1)+y−p) (10)

where , are parameters and is fixed. Since and have the same effect as bifurcation parameters we set . We consider several fixed y-values and the effect of varying (cf. Section 3.2) in each case. There are either one, two or three equilibrium points for (2.2). Equilibrium points satisfy and lie on the critical manifold, i.e. we have to solve

 0=x1(x1−1)(110−x1)+¯p (11)

for . We find that there are three equilibria for approximately , two equilibria on the boundary of this interval and one equilibrium otherwise. The Jacobian of (2.2) at an equilibrium is

 A(x1)=(01150(1−22x1+30x21)s5)

Direct calculation yields that for the single equilibrium is a saddle. 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 manifold , and . For the matrix is singular of rank 1 if and only if which occurs for the fold points . Hence the equilibria 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 values and . Note that by symmetry we can reduce to studying a single fold point. In the limit (corresponding to the case of a “standing wave”) the saddle-node bifurcation point becomes more degenerate with nilpotent.

Our next goal is to investigate global bifurcations of (2.2); we start with homoclinic orbits. For it is easy to see that (2.2) is a Hamiltonian system:

 x′1 = ∂H∂x2=x2 x′2 = −∂H∂x1=15(−x1(x1−1)(110−x1)−¯p) (12)

with Hamiltonian function

 H(x1,x2)=12x22−(x1)2100+11(x1)3150−(x1)420+x1¯p5 (13)

We will use this Hamiltonian formulation later on to describe the geometry of homoclinic orbits for slow wave speeds. Assume that is chosen so that (2.2) has a homoclinic orbit . We are interested in perturbations with and note that in this case the divergence of (2.2) is . Hence the vector field is area expanding everywhere. The homoclinic orbit breaks for and no periodic orbits are created. Note that this scenario does not apply to the full three-dimensional system as the equilibrium 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 -parameter space that represent heteroclinic connections of the fast subsystem. The main motivation is the decomposition of trajectories in the full system into slow and fast segments. Concatenating fast heteroclinic segments and slow flow segments can yield homoclinic orbits of the full system [28, 4, 32, 38]. We describe a numerical strategy to detect heteroclinic connections in the fast subsystem and continue them in parameter space. Suppose that so that (2.2) has three hyperbolic equilibrium points , and . We denote by the unstable and by the stable manifold of . The same notation is also used for and tangent spaces to and are denoted by and . Recall that is a source and shall not be of interest to us for now. Define the cross section by

 Σ={(x1,x2)∈R2:x1=xl+xr2}.

We use forward integration of initial conditions in and backward integration of initial conditions in to obtain trajectories and respectively. We calculate their intersection with and define

 γl(¯p,s):=γ+∩Σ,γr(¯p,s):=γ−∩Σ

We compute the functions and for different parameter values of numerically. Heteroclinic connections occur at zeros of the function

 h(¯p,s):=γl(¯p,s)−γr(¯p,s)

Once we find a parameter pair such that , these parameters can be continued along a curve of heteroclinic connections in parameter space by solving the root-finding problem for either or fixed and small. We use this method later for different fixed values of to compute heteroclinic connections in the fast subsystem in 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 and and approximately form a “”. From Figure 3 we conjecture that there exists a double heteroclinic orbit for .

Remarks: If we fix our initial change of variable becomes 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 heteroclinic connections of the fast subsystem can be used to prove the existence of homoclinic orbits to the unique saddle equilibrium (cf. [32]). Note that the existence of the heteroclinics in the fast subsystem was proved in a special case analytically [1] but Figure 3 is - to the best of our knowledge - the first explicit computation of where fast subsystem heteroclinics are located. The paper [36] develops a method for finding heteroclinic connections by the same basic approach we used, i.e. defining a codimension one hyperplane that separates equilibrium points.

Figure 3 suggests that there exists a double heteroclinic connection for . Observe that the Hamiltonian in our case is where the function is:

 V(x1)=px15−(x1)2100+11(x1)3150−(x1)420

The solution curves of (2.2) are given by . The structure of the solution curves entails symmetry under reflection about the -axis. Suppose and recall that we denoted the two saddle points of (2.2) by and and that their location depends on . Therefore, we conclude that the two saddles and must have a heteroclinic connection if they lie on the same energy level, i.e. they satisfy . This equation can be solved numerically to very high accuracy.

Proposition 1.

The fast subsystem of the FitzHugh-Nagumo equation for has a double heteroclinic connection for . Given a particular value there exists a double heteroclinic connection for in the fast subsystem lying in the plane .

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 -coordinate is a faster variable than . In particular, the small periodic orbits generated in the Hopf bifurcation lie almost in the plane . Hence to analyze this region we set to transform the FitzHugh-Nagumo equation (1.2) into a system with 2 slow and 1 fast variable:

 ˙x1 = ¯x2 ϵ2˙¯x2 = 15(sϵ¯x2−x1(x1−1)(110−x1)+y−p) (14) ˙y = 1s(x1−y)

Note that (2.3) corresponds to the FitzHugh-Nagumo equation in the form (cf. (4)):

 {uτ=5ϵ2uxx+f(u)−w+pwτ=ϵ(u−w) (15)

Therefore the transformation can be viewed as a rescaling of the diffusion strength by . We introduce a new independent small parameter and then let . This assumes that terms do not vanish in this limit, yielding the diffusion free system. Then the slow manifold of (2.3) is:

 S0={(x1,¯x2,y)∈R3:¯x2=1sϵ(f(x1)−y+p)} (16)
Proposition 2.

Following time rescaling by , the slow flow of system (2.3) on in the variables is given by

 ϵ˙x1 = f(x1)−y+p ˙y = x1−y (17)

In the variables the vector field (2) becomes

 ˙x1 = ¯x2 ϵ˙¯x2 = −1s2(x1−f(x1)−p)+¯x2s(f′(x1)−ϵ) (18)

Remark: The reduction to equations (2)-(2) suggests that (2.3) is a three time-scale system. Note however that (2.3) is not given in the three time-scale form for and . The time-scale separation in (2)-(2) results from the singular dependence of the critical manifold ; see (16).

Proof.

(of Proposition 2) Use the defining equation for the slow manifold (16) and substitute it into . A rescaling of time by under the assumption that yields the result (2). To derive (2) differentiate the defining equation of with respect to time:

 ˙¯x2=1sϵ(˙x1f′(x1)−˙y)=1sϵ(¯x2f′(x1)−˙y)

The equations and yield the equations (2). ∎

Before we start with the analysis of (2) we note that detailed bifurcation calculations for (2) exist. For example, Rocsoreanu et al. [44] give a detailed overview on the FitzHugh equation (2) and collect many relevant references. Therefore we shall only state the relevant bifurcation results and focus on the fast-slow structure and canards. Equation (2) has a critical manifold given by which coincides with the critical manifold of the full FitzHugh-Nagumo system (1.2). Formally it is located in but we still denote it by . Recall that the fold points are located at

 x1,±=130(11±√91)or numerically:x1,+≈0.6846,x1,−≈0.0487

Also recall that the y-nullcline passes through the fold points at:

 p−≈0.0511andp+≈0.5584

We easily find that supercritical Hopf bifurcations are located at the values

 pH,±(ϵ)=20576750±√11728171182250000−359ϵ1350+509ϵ22700−ϵ327 (19)

For the case we get and . The periodic orbits generated in the Hopf bifurcations exist for . Observe also that ; so the Hopf bifurcations of (2) coincide in the singular limit with the fold bifurcations in the one-dimensional slow flow (8). We are also interested in canards in the system and calculate a first order asymptotic expansion for the location of the maximal canard in (2) following [37]; recall that trajectories lying in the intersection of attracting and repelling slow manifolds are called maximal canards. We restrict to canards near the fold point .

Proposition 3.

Near the fold point the maximal canard in parameter space is given by:

 p(ϵ)=x1,−−c(x1,−)+58ϵ+O(ϵ3/2)
Proof.

Let and consider the shifts

 x1→x1+x1,−,¯y→¯y+c(x1,−),p→p+x1,−−c(x1,−)

to translate the equilibrium of (2) to the origin when . This gives

 x′1 = x21(√9110−x1)−¯y=¯f(x1,¯y) y′ = ϵ(x1−¯y−p)=ϵ(¯g(x1,¯y)−p) (20)

Now apply Theorem 3.1 in [37] to find that the maximal canard of (2.3) is given by:

 p(ϵ)=58ϵ+O(ϵ3/2)

Shifting the parameter back to the original coordinates yields the result. ∎

If we substitute in the previous asymptotic result and neglect terms of order then the maximal canard is predicted to occur for which is right after the first supercritical Hopf bifurcation at . Therefore we expect that there exist canard orbits evolving along the middle branch of the critical manifold in the full FitzHugh-Nagumo equation. Maximal canards are part of a process generally referred to as canard explosion [10, 39, 7]. In this situation the small periodic orbits generated in the Hopf bifurcation at undergo a transition to relaxation oscillations within a very small interval in parameter space. A variational integral determines whether the canards are stable [39, 26].

Proposition 4.

The canard cycles generated near the maximal canard point in parameter space for equation (2) are stable.

Proof.

Consider the differential equation (2) in its transformed form (2.3). Obviously this will not affect the stability analysis of any limit cycles. Let and denote the two smallest -coordinates of the intersection between

 ¯C0:={(x1,¯y)∈R2:¯y=√9110x21−x31=ϕ(x1)}

and the line . Geometrically represents a point on the left branch and a point on the middle branch of the critical manifold . Theorem 3.4 in [39] tells us that the canards are stable cycles if the function

 R(h)=∫xm(h)xl(h)∂¯f∂x1(x1,ϕ(x1))ϕ′(x1)¯g(x1,ϕ(x1))dx1

is negative for all values where is the second fold point of besides . In our case we have

 R(h)=∫xm(h)xl(h)(√915x1−3x21)2x−√9110x21+x31dx

with and . Figure 4 shows a numerical plot of the function for the relevant values of which confirms the required result.

Remark: We have computed an explicit algebraic expression for with a computer algebra system. This expression yields for , confirming that is decreasing.

As long as we stay on the critical manifold of the full system, the analysis of the bifurcations and geometry of (2) give good approximations to the dynamics of the FitzHugh-Nagumo equation because the rescaling leaves the plane invariant. Next we use the dynamics of the -coordinate in system (2) to obtain better insight into the dynamics when . The critical manifold of (2) is:

 D0={(x1,¯x2)∈R2:s¯x2c′(x1)=x1−c(x1)}

We are interested in the geometry of the periodic orbits shown in Figure 5 that emerge from the Hopf bifurcation at . Observe that the amplitude of the orbits in the direction is much larger that than in the -direction. Therefore we predict only a single small excursion in the direction for slightly larger than as shown in Figures 5(a) and 5(c). The wave speed changes the amplitude of this excursion with a smaller wave speed implying a larger excursion. Hence equation (2) is expected to be a very good approximation for periodic orbits in the FitzHugh-Nagumo equation with fast wave speeds. Furthermore the periodic orbits show two excursions in the relaxation regime after the canard explosion; see Figure 5(b).

3 The Full System

3.1 Hopf Bifurcation

The characteristic polynomial of the linearization of the FitzHugh-Nagumo equation (1.2) at its unique equilibrium point is

 P(λ)=ϵ5s+(−ϵs−λ)(−150+11x∗125−3(x∗1)25−sλ5+λ2)

Denoting , a necessary condition for to have pure imaginary roots is that . The solutions of this equation can be expressed parametrically as a curve :

 s(x∗1)2 = 50ϵ(ϵ−1)1+10ϵ−22x∗1+30(x∗1)2 p(x∗1) = (x∗1)3−1.1(x∗1)2+1.1 (21)
Proposition 5.

In the singular limit the U-shaped bifurcation curves of the FitzHugh-Nagumo equation have vertical asymptotes given by the points and and a horizontal asymptote given by . Note that at the equilibrium point passes through the two fold points.

Proof.

The expression for in (3.1) is positive when . For values of between the roots of , in (3.1) as . The values of and in the proposition are approximations to the value of in (3.1) at the roots of . As , solutions of the equation in (3.1) yield values of that tend to one of the two roots of . The result follows. ∎

The analysis of the slow subsystems (2) and (2) gives a conjecture about the shape of the periodic orbits in the FitzHugh-Nagumo equation. Consider the parameter regime close to a Hopf bifurcation point. From (2) we expect one part of the small periodic orbits generated in the Hopf bifurcation to lie close to the slow manifolds and . Using the results about equation (2) we anticipate the second part to consist of an excursion in the direction whose length is governed by the wave speed . Figure 6 shows a numerical continuation in MatCont [22] of the periodic orbits generated in a Hopf bifurcation and confirms the singular limit analysis for small amplitude orbits.

Furthermore we observe from comparison of the and coordinates of the periodic orbits in Figure 6(b) that orbits tend to lie close to the plane defined by . More precisely, the diameter of the periodic orbits is observed to be in this case. This indicates that the rescaling of Section 2.3 can help to describe the system close to the U-shaped Hopf curve. Note that it is difficult to check whether this observation of an -diameter in the -coordinate persists for values of since numerical continuation of canard-type periodic orbits is difficult to use for smaller .

In contrast to this, it is easily possible to compute the U-shaped Hopf curve using numerical continuation for very small values of . We have used this possibility to track two generalized Hopf bifurcation points in three parameters . The U-shaped Hopf curve has been computed by numerical continuation for a mesh of parameter values for between and using MatCont [22]. The two generalized Hopf points on the left half of the U-curve were detected as codimension two points during each continuation run. The results of this “three-parameter continuation” are shown in Figure 7.

The two generalized Hopf points depend on and we find that their singular limits in -parameter space are approximately:

 GH01≈(p=0.171,s=0)andGH02≈(p=0.051,s=3.927)

We have not found a way to recover these special points from the fast-slow decomposition of the system. This suggests that codimension two bifurcations are generally diffcult to recover from the singular limit of fast-slow systems.

Furthermore the Hopf bifurcations for the full system on the left half of the U-curve are subcritical between and and supercritical otherwise. For the transformed system (2.3) with two slow and one fast variable we observed that in the singular limit (2) for the Hopf bifurcation is supercritical. In the case of the periodic orbits for (1.2) and (2) exist in overlapping regions for the parameter between the -values of and . This result indicates that (2.3) can be used to describe periodic orbits that will interact with the homoclinic C-curve.

3.2 Homoclinic Orbits

In the following discussion we refer to “the” C-shaped curve of homoclinic bifurcations of system (1.2) as the parameters yielding a “single-pulse” homoclinic orbit. The literature as described in Section 1.2 shows that close to single-pulse homoclinic orbits we can expect multi-pulse homoclinic orbits that return close to the equilibrium point multiple times. Since the separation of slow manifolds is exponentially small, homoclinic orbits of different types will always occur in exponentially thin bundles in parameter space. Values of are small enough that the parameter region containing all the homoclinic orbits will be indistinguishable numerically from “the” C-curve that we locate.

The history of proofs of the existence of homoclinic orbits in the FitzHugh-Nagumo equation is quite extensive. The main step in their construction is the existence of a “singular” homoclinic orbit . We consider the case when the fast subsystem has three equilibrium points which we denote by , and . Recall that coincides with the unique equilibrium of the full system for . A singular homoclinic orbit is always constructed by first following the unstable manifold of in the fast subsystem given by .

First assume that . In this case the Hamiltonian structure - see Section 2.2 and equation (2.2) - can be used to show the existence of a singular homoclinic orbit. Figure 8 shows level curves for various values of . The double heteroclinic connection can be calculated directly using Proposition 1 and solving for .

Proposition 6.

There exists a singular double heteroclinic connection in the FitzHugh-Nagumo equation for and .

Techniques developed in [45] show that the singular homoclinic orbits existing for and must persist for perturbations of small positive wave speed and sufficiently small . These orbits are associated to the lower branch of the C-curve. The expected geometry of the orbits is indicated by their shape in the singular limit shown in Figure 8. The double heteroclinic connection is the boundary case between the upper and lower half of the C-curve. It remains to analyze the singular limit for the upper half. In this case, a singular homoclinic orbit is again formed by following the unstable manifold of when it coincides with the equilibrium but now we check whether it forms a heteroclinic orbit with the stable manifold of . Then we follow the slow flow on and return on a heteroclinic connection to for a different y-coordinate with and . From there we connect back via the slow flow. Using the numerical method described in Section 2.2 we first set ; note that the location of depends on the value of the parameter . The task is to check when the system

 x′1 = x2 x′2 = 15(f(x1)+y−p) (22)

has heteroclinic orbits from to with . The result of this computation is shown in Figure 9 as the red curve. We have truncated the result at . In fact, the curve in Figure 9 can be extended to . Obviously we should view this curve as an approximation to the upper part of the C-curve.

If the connection from back to occurs with vertical coordinate , it is a trajectory of system (3.2) with . Figure 9 shows values of at which these heteroclinic orbits exist for . An intersection between a red and a blue curve indicates a singular homoclinic orbit. Further computations show that increasing the value of slowly beyond yields intersections everywhere along the red curve in Figure 9. Thus the values of on the homoclinic orbits are expected to grow as increases along the upper branch of the C-curve. Since there cannot be any singular homoclinic orbits for we have to find the intersection of the red curve in Figure 9 with the vertical line . Using the numerical method to detect heteroclinic connections gives:

Proposition 7.

The singular homoclinic curve for positive wave speed terminates at and on the right and at and on the left.

In -parameter space define the points:

 A=(p∗,0),B=(p−,0),C=(p−,s∗) (23)

In Figure 10 we have computed the homoclinic C-curve for values of between and . Together with the singular limit analysis above, this yields strong numerical evidence for the following conjecture:

Conjecture 1.

The C-shaped homoclinic bifurcation curves converge to the union of the segments and as .

Remark 1: Figure 4 of Krupa, Sandstede and Szmolyan [38] shows a “wedge” that resembles shown in Figure 10. The system that they study sets and varies with . For and , the equilibrium point is located at the origin and the fast subsystem with has a double heteroclinic connection at to the saddle equilibrium . The techniques developed in [38] use this double heteroclinic connection as a starting point. Generalizations of the results in [38] might provide a strategy to prove Conjecture 1 rigorously, a possibility that we have not yet considered. However, we think that 1-homoclinic orbits in the regime we study come in pairs and that the surface of 1-homoclinic orbits in space differs qualitatively from that described by Krupa, Sandstede and Szmolyan.

Remark 2: We have investigated the termination or turning mechanism of the C-curve at its upper end. The termination points shown in Figure 1 have been obtained by a different geometric method. It relies on the observation that, in addition to the two fast heteroclinic connections, we have to connect near back to the equilibrium point to form a homoclinic orbit; the two heteroclinic connections might persist as intersections of suitable invariant manifolds but we also have to investigate how the flow near interacts with the stable manifold . These results will be reported elsewhere, but we note here that .

The numerical calculations of the C-curves for are new. Numerical continuation using the boundary value methods implemented in AUTO [14] or MatCont [22] becomes very difficult for these small values of [5]. Even computing with values using boundary value methods is a numerically challenging problem. The method we have used does not compute the homoclinic orbits themselves while it locates the homoclinic C-curve accurately in parameter space. To motivate our approach consider Figure 11 which shows the unstable manifold for different values of and fixed . We observe that homoclinic orbits can only exist near two different wave speeds and which define the parameters where or . Figure 11 displays how changes as varies for the fixed value . If differs from the points and that define the lower and upper branches of the C-curve for the given value of , then increases rapidly on away from . The changes in sign of on identify values of with homoclinic orbits. The two splitting points that mark these sign changes are visible in Figure 11. Since trajectories close to the slow manifolds separate exponentially away from them, we are able to assess the sign of unambiguously on trajectories close to the slow manifold and find small intervals and that contain the values of for which there are homoclinic orbits.

The geometry of the orbits along the upper branch of the C-curve is obtained by approximating it with two fast singular heteroclinic connections and parts of the slow manifolds and ; this process has been described several times in the literature when different methods were used to prove the existence of “fast waves” (see e.g. [29, 4, 32]).

4 Conclusions

Our results are summarized in the singular bifurcation diagram shown in Figure