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

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 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.


Fast-Slow Systems
Fast-slow systems of ordinary differential equations (ODEs) have the general form: where x ∈ R m , y ∈ R n 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 algebraic constraint f = 0 defines the critical manifold C 0 = {(x, y) ∈ R m × R n : f (x, y, 0) = 0}.
Where D x f (p) is nonsingular, the implicit function theorem implies that there exists a map h(x) = y parametrizing C 0 as a graph. This yields the implicitly defined vector fieldẏ = g(h(y), y, 0) on C 0 called the slow flow.
We can change (1) to the fast time scale t = τ /ǫ and let ǫ → 0 to obtain the second possible singular limit system We call the vector field (2) parametrized by the slow variables y 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 C 0 normally hyperbolic there exists a family of invariant manifolds C ǫ for the flow (1). The manifolds are at a distance O(ǫ) from C 0 and the flows on them converge to the slow flow on C 0 as ǫ → 0. Points p ∈ C 0 where D x f (p) is singular are referred to as fold points 1 .
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.

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: 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: where f a (u) = u(u − a)(1 − u) and p, γ, δ and a 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 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 xx = u ′′ and w τ = sw ′ . Set v = u ′ and substitute into (4) to obtain the system: 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 analyzed using several different methods. For example, it has been proved that (5) 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 (5) 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 (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 codimension 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 employed 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: We shall fix those values to allow comparison of our results with theirs. Hence we also write f 1/10 (u) = f (u). Changing from the fast time t to the slow time τ and relabeling variables x 1 = u, x 2 = v and y = w we get: From now on we refer to (6) as "the" FitzHugh-Nagumo equation. Investigating bifurcations in the (p, s) 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 (6) has the symmetry 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 ǫ ≤ 10 −3 . 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.   (6). Hopf bifurcations are shown in green, saddle-node of limit cycles (SNLC) are shown in blue and GH indicates a generalized 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 two black dots (see Section 3.2). Only part of the parameter space is shown because of the symmetry (7). The homoclinic curve has been thickened to indicate that multipulse homoclinic orbits exist very close to single pulse ones (see [15]).
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 ǫ → 0.

Remark:
We have also investigated the termination mechanism of the Cshaped 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.

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: Therefore C 0 is a cubic curve in the coordinate plane x 2 = 0. The parameter p moves the cubic up and down inside this plane. The critical points of the cubic are solutions of c ′ (x 1 ) = 0 and are given by: The points x 1,± are fold points with |c ′′ (x 1,± )| = 0 since C 0 is a cubic polynomial with distinct critical points. The fold points divide C 0 into three segments We denote the associated slow manifolds by C l,ǫ , C m,ǫ and C r,ǫ . There are two possibilities to obtain the slow flow. One way is to solve c(x 1 ) = y for x 1 and substitute the result into the equationẏ = 1 s (x 1 − y). Alternatively differentiating y = c(x 1 ) implicitly with respect to τ yieldsẏ =ẋ 1 c ′ (x 1 ) and therefore One can view this as a projection of the slow flow, which is constrained to the critical manifold in R 3 , onto the x 1 -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 y = x 1 only in a single point x * 1 which is the unique equilibrium of the slow flow (8). Observe that q = (x * 1 , 0, x * 1 ) is also the unique 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 0 is to look at the sign of (x 1 − y). The situation is illustrated in Figure 2. where the subscripts indicate the fold point at which each equilibrium is located.
The singular time-rescalingτ = sc ′ (x 1 )/τ of the slow flow yields the desingularized slow flow Time is reversed by this rescaling on C l and C r since s > 0 and c ′ (x 1 ) is negative on these branches. The desingularized slow flow (9) is smooth and has no bifurcations as p is varied.

The Fast Subsystem
The key component of the fast-slow analysis for the FitzHugh-Nagumo equation is the two-dimensional fast subsystem 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 2 = 0 and lie on the critical manifold, i.e. we have to solve for x 1 . We find that there are three equilibria for approximatelyp l = −0.1262 < p < 0.0024 =p r , two equilibria on the boundary of this p interval and one equilibrium otherwise. The Jacobian of (10) at an equilibrium is Direct calculation yields that for p ∈ [p l ,p r ] 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 C l , C m and C r . For s > 0 the matrix A is singular of rank 1 if and only if 30x 2 1 − 22x 1 + 1 = 0 which occurs for the fold points x 1,± . 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 valuesp l andp r . Note that by symmetry we can reduce to studying a single fold 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 1 ) nilpotent.
Our next goal is to investigate global bifurcations of (10); we start with homoclinic orbits. For s = 0 it is easy to see that (10) is a Hamiltonian system: with Hamiltonian function We will use this Hamiltonian formulation later on to describe the geometry of homoclinic orbits for slow wave speeds. Assume thatp is chosen so that (12) has a homoclinic orbit x 0 (t). We are interested in perturbations with s > 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 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 thatp ∈ (p l ,p r ) so that (10) has three hyperbolic equilibrium points x l , x m and x r . We denote by W u (x l ) the unstable and by W s (x l ) the stable manifold of x l . The same notation is also used for x r and tangent spaces to W s (.) and W u (.) are denoted by T s (.) and T u (.). Recall that x m is a source and shall not be of interest to us for now. Define the cross section Σ by We use forward integration of initial conditions in T u (x l ) and backward integration of initial conditions in T s (x r ) to obtain trajectories γ + and γ − respectively. We calculate their intersection with Σ and define We compute the functions γ l and γ r for different parameter values of (p, s) numerically. Heteroclinic connections occur at zeros of the function Once we find a parameter pair (p 0 , s 0 ) such that h(p 0 , s 0 ) = 0, these parameters can be continued along a curve of heteroclinic connections in (p, s) parameter space by solving the root-finding problem h(p 0 + δ 1 , s 0 + δ 2 ) = 0 for either δ 1 or δ 2 fixed and small. We use this method later for different fixed values of y 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 top l andp r and approximately form a "V ". From Figure 3 we conjecture that there exists a double heteroclinic orbit forp ≈ −0.0622. 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 heteroclinic 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 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 H that separates equilibrium points. Figure 3 suggests that there exists a double heteroclinic connection for s = 0. Observe that the Hamiltonian in our case is H( The solution curves of (12) are given by x 2 = ± 2(const. − V (x 1 )). The structure of the solution curves entails symmetry under reflection about the x 1axis. Supposep ∈ [p l ,p r ] and recall that we denoted the two saddle points of (10) by x l and x r and that their location depends onp. Therefore, we conclude that the two saddles x l and x r must have a heteroclinic connection if they lie on the same energy level, i.e. they satisfy V (x l ) − V (x r ) = 0. This equation 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 forp =p * ≈ −0.0619259. Given a particular value y = y 0 there exists a double heteroclinic connection for p = p * + y 0 in the fast subsystem lying in the plane y = y 0 .

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 2 -coordinate is a faster variable than x 1 . In particular, the small periodic orbits generated in the Hopf bifurcation lie almost in the plane x 2 = 0. Hence to analyze this region we setx 2 = x 2 /ǫ to transform the FitzHugh-Nagumo equation (6) into a system with 2 slow and 1 fast variable: Note that (14) corresponds to the FitzHugh-Nagumo equation in the form (cf. (4)): Therefore the transformationx 2 = x 2 /ǫ can be viewed as a rescaling of the 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 0 of (14) is: Proposition 2. Following time rescaling by s, the slow flow of system (14) on S 0 in the variables (x 1 , y) is given by In the variables (x 1 ,x 2 ) the vector field (17) becomeṡ Remark: The reduction to equations (17)- (18) suggests that (14) is a three time-scale system. Note however that (14) is not given in the three time-scale form (ǫ 2ż 1 , ǫż 2 ,ż 3 ) = (h 1 (z), h 2 (z), h 3 (z)) for z = (z 1 , z 2 , z 3 ) ∈ R 3 and h i : 1, 2, 3). The time-scale separation in (17)-(18) results from the singular 1/ǫ dependence of the critical manifold S 0 ; see (16).
Proof. (of Proposition 2) Use the defining equation for the slow manifold (16) and substitute it intoẋ 1 =x 2 . A rescaling of time by t → st under the assumption that s > 0 yields the result (17). To derive (18) differentiate the defining equation of S 0 with respect to time: The equationsẏ = 1 s (x 1 − y) and y = −sǫx 2 + f (x 1 ) + p yield the equations (18).
Before we start with the analysis of (17) we note that detailed bifurcation calculations for (17) exist. For example, Rocsoreanu et al. [44] give a detailed overview on the FitzHugh equation (17) and collect many relevant references. Therefore we shall only state the relevant bifurcation results and focus on the fast-slow structure and canards. Equation (17) has a critical manifold given by y = f (x 1 ) + p = c(x 1 ) which coincides with the critical manifold of the full FitzHugh-Nagumo system (6). Formally it is located in R 2 but we still denote it by C 0 . Recall that the fold points are located at Also recall that the y-nullcline passes through the fold points at: p − ≈ 0.0511 and p + ≈ 0.5584 We easily find that supercritical Hopf bifurcations are located at the values For the case ǫ = 0.01 we get p H,− (0.01) ≈ 0.05632 and p H,+ (0.01) ≈ 0.55316. The periodic orbits generated in the Hopf bifurcations exist for p ∈ (p H,− , p H,+ ). Observe also that p H,± (0) = p ± ; so the Hopf bifurcations of (17) 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 (17) 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 (x 1,− , c(x 1,− )).
Proposition 3. Near the fold point (x 1,− , c(x 1,− )) the maximal canard in (p, ǫ) parameter space is given by: Proof. Letȳ = y − p and consider the shifts to translate the equilibrium of (17) to the origin when p = 0. This gives Now apply Theorem 3.1 in [37] to find that the maximal canard of (20) is given by: Shifting the parameter p back to the original coordinates yields the result.
If we substitute ǫ = 0.01 in the previous asymptotic result and neglect terms of order O(ǫ 3/2 ) then the maximal canard is predicted to occur for p ≈ 0.05731 which is right after the first supercritical Hopf bifurcation at p H,− ≈ 0.05632. Therefore we expect that there exist canard orbits evolving along the middle branch of the critical manifold C m,0.01 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 p = p H,− 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]. Proof. Consider the differential equation (17) in its transformed form (20). Obviously this will not affect the stability analysis of any limit cycles. Let x l (h) and x m (h) denote the two smallest x 1 -coordinates of the intersection between 15 is the second fold point ofC 0 besides x 1 = 0. In our case we have
As long as we stay on the critical manifold C 0 of the full system, the analysis of the bifurcations and geometry of (17) give good approximations to the dynamics of the FitzHugh-Nagumo equation because the rescaling x 2 = ǫx 2 leaves the plane x 2 = 0 invariant. Next we use the dynamics of thex 2 -coordinate in system (18) to obtain better insight into the dynamics when x 2 = 0. The critical manifold D 0 of (18) is: We are interested in the geometry of the periodic orbits shown in Figure 5 that emerge from the Hopf bifurcation at p H,− . Observe that the amplitude of the orbits in the x 1 direction is much larger that than in the x 2 -direction. Therefore we predict only a single small excursion in the x 2 direction for p slightly larger than p H,− as shown in Figures 5(a) and 5(c). The wave speed changes the amplitude of this x 2 excursion with a smaller wave speed implying a larger excursion. Hence equation (17) 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 x 2 excursions in the relaxation regime after the canard explosion; see Figure 5(b).  (18). Note that here x 2 = ǫx 2 is shown. Orbits have been obtained by direct forward integration for ǫ = 0.01.

Hopf Bifurcation
The characteristic polynomial of the linearization of the FitzHugh-Nagumo equation (6) at its unique equilibrium point is Denoting P (λ) = c 0 + c 1 λ + c 2 λ 2 + c 3 λ 3 , a necessary condition for P to have pure imaginary roots is that c 0 = c 1 c 2 . The solutions of this equation can be expressed parametrically as a curve (p(x * 1 ), s(x * 1 )): The analysis of the slow subsystems (17) and (18) 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 (17) we expect one part of the small periodic orbits generated in the Hopf bifurcation to lie close to the slow manifolds C l,ǫ and C m,ǫ . Using the results about equation (18) we anticipate the second part to consist of an excursion in the x 2 direction whose length is governed by the wave speed s. 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 x 1 and x 2 coordinates of the periodic orbits in Figure 6(b) that orbits tend to lie close to the plane defined by x 2 = 0. More precisely, the x 2 diameter of the periodic orbits is observed to be O(ǫ) 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 O(ǫ)-diameter in the x 2 -coordinate persists for values of ǫ < 0.01 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 (p, s, ǫ). The U-shaped Hopf curve has been computed by numerical continuation for a mesh of parameter values for ǫ between 10 −2 and 10 −7 using MatCont [22]. The two generalized Hopf points GH ǫ 1,2 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 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 GH ǫ 1 and GH ǫ 2 and supercritical otherwise. For the transformed system (14) with two slow and one fast variable we observed that in the singular limit (17) for ǫ 2 → 0 the Hopf bifurcation is supercritical. In the case of ǫ = 0.01 the periodic orbits for (6) and (17)  . This result indicates that (14) can be used to describe periodic orbits that will interact with the homoclinic C-curve.

Homoclinic Orbits
In the following discussion we refer to "the" C-shaped curve of homoclinic bifurcations of system (5) 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 C ·,ǫ is exponentially small, homoclinic orbits of different types will always occur in exponentially thin bundles in parameter space. Values of ǫ < 0.005 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 γ 0 . We consider the case when the fast subsystem has three equilibrium points which we denote by x l ∈ C l , x m ∈ C m and x r ∈ C r . Recall that x l coincides with the unique equilibrium q = (x * 1 , 0, x * 1 ) of the full system for p < p − . A singular homoclinic orbit is always constructed by first following the unstable manifold of x l in the fast subsystem given by y = x * 1 .  First assume that s = 0. In this case the Hamiltonian structure -see Section 2.2 and equation (12) -can be used to show the existence of a singular homoclinic orbit. Figure 8 shows level curves H(x 1 , x 2 ) = H(x * 1 , 0) for various values of p. The double heteroclinic connection can be calculated directly using Proposition 1 and solving x * 1 +p * = p for p. Proposition 6. There exists a singular double heteroclinic connection in the FitzHugh-Nagumo equation for s = 0 and p ≈ −0.246016 = p * .
Techniques developed in [45] show that the singular homoclinic orbits existing for s = 0 and p ∈ (p * , p − ) 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 Ccurve. 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 x l when it coincides with the equilibrium q = (x * 1 , 0, x * 1 ) but now we check whether it forms a heteroclinic orbit with the stable manifold of x r . Then we follow the slow flow on C r and return on a heteroclinic connection to C l for a different y-coordinate with y > x * 1 and y < c(x 1,+ ) = f (x 1 ) + p. From there we connect back via the slow flow. Using the numerical method described in Section 2.2 we first set y = x * 1 ; note that the location of q depends on the value of the parameter p. The task is to check when the system has heteroclinic orbits from C l to C r with y = x * 1 . The result of this computation is shown in Figure 9 as the red curve. We have truncated the result at p = −0.01. In fact, the curve in Figure 9 can be extended to p = p * . Obviously we should view this curve as an approximation to the upper part of the C-curve. If the connection from C r back to C l occurs with vertical coordinate x * 1 +v, it is a trajectory of system (22) with y = x * 1 + v. Figure 9 shows values of (p, s) at which these heteroclinic orbits exist for v = 0.125, 0.12, 0.115. An intersection between a red and a blue curve indicates a singular homoclinic orbit. Further computations show that increasing the value of v slowly beyond 0.125 yields intersections everywhere along the red curve in Figure 9. Thus the values of v on the homoclinic orbits are expected to grow as s increases along the upper branch of the C-curve. Since there cannot be any singular homoclinic orbits for p ∈ (p − , p + ) we have to find the intersection of the red curve in Figure 9 with the vertical line p = p − . Using the numerical method to detect heteroclinic connections gives: Proposition 7. The singular homoclinic curve for positive wave speed terminates at p = p − and s ≈ 1.50815 = s * on the right and at p = p * and s = 0 on the left.
In (p, s)-parameter space define the points: In Figure 10 we have computed the homoclinic C-curve for values of ǫ between 10 −2 and 5 · 10 −5 . 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 AB and AC as ǫ → 0.
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 p = 0 and varies a with a ≈ 1/2. For a = 1/2 and p = 0, the equilibrium point q is located at the origin and the fast subsystem with y = 0 has a double heteroclinic connection at q to the saddle equilibrium (1, 0, 0) ∈ C r . 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 (p, s, ǫ) 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 C l back to the equilibrium point q 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 C l,ǫ interacts with the stable manifold W s (q). These results will be reported elsewhere, but we note here that p turn (ǫ) → p − .
The numerical calculations of the C-curves for ǫ ≤ 10 −3 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 ǫ = O(10 −2 ) 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 W u (q) for different values of s and fixed p. We observe that homoclinic orbits can only exist near two different wave speeds s 1 and s 2 which define the parameters where W u (q) ⊂ W s (C l,ǫ ) or W u (q) ⊂ W s (C r,ǫ ).  11 displays how W u (q) changes as s varies for the fixed value p = 0.05. If s differs from the points s 1 and s 2 that define the lower and upper branches of the C-curve for the given value of p, then |x 1 | increases rapidly on W u (q) away from q. The changes in sign of x 1 on W u (q) identify values of s 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 x 1 unambiguously on trajectories close to the slow manifold and find small intervals (p, s 1 ± 10 −15 ) and (p, s 2 ± 10 −15 ) that contain the values of s 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 C r,ǫ and C l,ǫ ; 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]).

Conclusions
Our results are summarized in the singular bifurcation diagram shown in Figure 12. This figure shows information obtained by a combination of fast-slow decompositions, classical dynamical systems techniques and robust numerical algorithms that work for very small values of ǫ. It recovers and extends to smaller values of ǫ the CU-structure described in [5] for the FitzHugh-Nagumo equation. The U-shaped Hopf curve was computed with an explicit formula, and the homoclinic C-curve was determined by locating transitions between different dynamical behaviors separated by the homoclinic orbits. All the results shown as solid lines in Figure 12 have been obtained by considering a singular limit. The lines AB and AC as well as the slow flow bifurcation follow from the singular limit ǫ → 0 yielding the fast and slow subsystems of the FitzHugh-Nagumo equation (6). The analysis of canards and periodic orbits have been obtained from equations (17) and (18) where the singular limit ǫ 2 → 0 was used (see Section 2.3). We have also shown the C-and U-curves in Figure 12 as dotted lines to orient the reader how the results from Proposition 5 and Conjecture 1 fit in.
We also observed that several dynamical phenomena are difficult to recover from the singular limit fast-slow decomposition. In particular, the codimension two generalized Hopf bifurcation does not seem to be observable from the singular limit analysis. Furthermore the homoclinic orbits can be constructed from the singular limits but it cannot be determined directly from the fast and slow subsystems that they are of Shilnikov-type.
The type of analysis pursued here seems to be very useful for other multiple time-scale problems involving multi-parameter bifurcation problems. In future work, we shall give a geometric analysis of the folding/turning mechanism of the homoclinic C-curve, a feature of this system we have not been able to determine directly from our singular limit analysis. That work relies upon new methods for calculating C l,ǫ and C r,ǫ which are invariant slow manifolds of "saddle-type" with both stable and unstable manifolds.
We end with brief historical remarks. The references cited in this paper  (17),(18) Figure 12: Sketch of the singular bifurcation diagram for the FitzHugh-Nagumo equation (6). The points A, B and C are defined in (23). The part of the diagram obtained from equations (17), (18) corresponds to the case "ǫ 2 = 0 and ǫ = 0 and small". In this scenario the canards to the right of p = p − are stable (see Proposition 4). The phase portrait in the upper right for equation (17) shows the geometry of a small periodic orbit generated in the Hopf bifurcation of (17). The two phase portraits below it show the geometry of these periodic orbits further away from the Hopf bifurcation for (17), (18). Excursions of the periodic orbits/canards for p > p − decrease for larger values of s. Note also that we have indicated as dotted lines the C-curve and the U-curve for positive ǫ to allow a qualitative comparison with Figure 1. discuss mathematical challenges posed by the FitzHugh-Nagumo equation, how these challenges have been analyzed and their relationship to general questions about multiple time-scale systems. Along the line AB in Figure 12 we encounter a perturbation problem regarding the persistence of homoclinic orbits that can be solved using Fenichel theory [45]. The point A marks the connection between fast and slow waves in (p, s)-parameter space which has been investigated in (ǫ, s)-parameter space in [38]. We view this codimension 2 connectivity as one of the key features of the FitzHugh-Nagumo system. The perturbation problem for homoclinic orbits close to the line AC was solved using several methods and was put into the context of multiple time-scale systems in [31,32], where the Exchange Lemma overcame difficulties in tracking W u (q) when it starts jumping away from C r,ǫ . This theory provides rigorous foundations that support our numerical computations and their interpretation.