Saddle Slow Manifolds and Canard Orbits in R4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{R}^{4}$\end{document} and Application to the Full Hodgkin–Huxley Model

Many physiological phenomena have the property that some variables evolve much faster than others. For example, neuron models typically involve observable differences in time scales. The Hodgkin–Huxley model is well known for explaining the ionic mechanism that generates the action potential in the squid giant axon. Rubin and Wechselberger (Biol. Cybern. 97:5–32, 2007) nondimensionalized this model and obtained a singularly perturbed system with two fast, two slow variables, and an explicit time-scale ratio ε. The dynamics of this system are complex and feature periodic orbits with a series of action potentials separated by small-amplitude oscillations (SAOs); also referred to as mixed-mode oscillations (MMOs). The slow dynamics of this system are organized by two-dimensional locally invariant manifolds called slow manifolds which can be either attracting or of saddle type. In this paper, we introduce a general approach for computing two-dimensional saddle slow manifolds and their stable and unstable fast manifolds. We also develop a technique for detecting and continuing associated canard orbits, which arise from the interaction between attracting and saddle slow manifolds, and provide a mechanism for the organization of SAOs in R4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{R}^{4}$\end{document}. We first test our approach with an extended four-dimensional normal form of a folded node. Our results demonstrate that our computations give reliable approximations of slow manifolds and canard orbits of this model. Our computational approach is then utilized to investigate the role of saddle slow manifolds and associated canard orbits of the full Hodgkin–Huxley model in organizing MMOs and determining the firing rates of action potentials. For ε sufficiently large, canard orbits are arranged in pairs of twin canard orbits with the same number of SAOs. We illustrate how twin canard orbits partition the attracting slow manifold into a number of ribbons that play the role of sectors of rotations. The upshot is that we are able to unravel the geometry of slow manifolds and associated canard orbits without the need to reduce the model.


Introduction
In a wide variety of models, including of semiconductor lasers [2,3], chemical reactions [4][5][6] and neurons [7][8][9], one finds large differences in time scales. For example, in neuron models, the action potential changes on a much faster time scale than ionic gating variables. This yields various types of complex oscillatory behavior such as spiking, bursting and mixed-mode oscillations (MMOs),which alternate between small-amplitude oscillations (SAOs) and large-amplitude oscillations (LAOs) [10]. Those phenomena can be modeled by a special class of dynamical systems called slow-fast systems, which have a group of fast and a group of slow variables, with a time-scale separation parameter ε. A common approach for studying slow-fast systems is geometric singular perturbation theory (GSPT) [11,12] which exploits the time-scale differences and approximates a slow-fast system by a composition of lower-dimensional slow and fast subsystems that are easier to analyze.
In the simplest case, slow-fast systems feature one fast and one slow variables. One famous example is the FitzHugh-Nagumo model [13,14], which is a simplified two-dimensional dynamical system that describes the evolution of the voltage potential and a single coupled slow gating variable. The voltage nullcline corresponds to the so-called critical manifold. The cubic shape of the critical manifold provides a mechanism for spike generation and relaxation oscillations. The notion relaxation oscillation was first introduced for the Van der Pol oscillator [15] to describe periodic orbits with slow and fast segments. Parameter variations may lead to transitions between an excitable rest state, small periodic orbits and large relaxation oscillations. The transformation of small periodic orbits of the FitzHugh-Nagumo model into large relaxation oscillations occurs within an exponentially small parameter range. In this parameter range, periodic orbits known as canard cycles follow the middle (repelling) branch of the critical manifold for a considerable amount of time. This phenomenon is called "canard explosion" and was first discovered by Benoît in 1981 [16]. The analysis of canards has since been prominent for understanding slow-fast systems.
Neuron models with one slow and two fast variables may feature high-frequency spikes separated by quiescent periods. This behavior is called "bursting". Rinzel [17] introduced a scheme for classifying bursting mechanisms based on the bifurcation structure of the so-called fast subsystem where the slow variable is treated as a parameter. Based on this structure, Izhikevich [18] attempted to classify all such bursting patterns. The critical manifold in bursting models with one slow and two fast variables is a one-dimensional manifold of equilibrium points of the fast subsystem, and it often takes a cubic shape. In this setting, canard orbits feature segments that stay close to the middle (saddle) branch of the critical manifold for a significant amount of time. Canard orbits play a key role in organizing bursting patterns; namely, they allow for an abrupt increase of amplitude in a very small parameter range. This creates a mechanism for spike adding [19] in bursting models [9,[20][21][22][23][24][25][26]. A different mechanism specifically for the transition between spiking and bursting has been related to so-called torus canards, which are orbit segments that connect the stable and unstable branches of periodic orbits of the fast subsystem [27][28][29][30][31].
Systems with one fast and two slow variables may generate a robust mechanism for MMOs. The critical manifolds of such systems are typically two-dimensional folded surfaces. Attracting and repelling sheets of a folded critical manifold meet at fold curves. According to well-established results from Fenichel [11,12], normally hyperbolic sheets of the critical manifold persist as two-dimensional slow manifolds, provided ε is sufficiently small. Intersections of attracting and repelling slow manifolds in this setting correspond to canard orbits. They provide a mechanism for SAOs in that canard orbits form boundaries between regions of different numbers of SAOs [32][33][34]. Global return mechanisms and canard-induced SAOs together generate MMOs for a wide range of mathematical and applied three-dimensional models [1,9,20,[35][36][37][38][39][40][41][42][43][44]. In support of theoretical work, numerical methods were developed to reveal the geometry of two-dimensional attracting and repelling slow manifolds, as well as canard orbits. These methods were employed to investigate numerous threedimensional examples [35,36,41,43,[45][46][47][48]. For a more detailed survey of MMOs in slow-fast systems see [10].
Mechanisms of bursting and MMOs in three-dimensional vector fields are relatively well understood. Nevertheless, there are many three-dimensional slow-fast systems [1,37,49,50] that are actually reductions of higher-dimensional models. For example, the Hodgkin-Huxley model is a four-dimensional vector field that was originally formed to explain the ionic mechanism that generates the action potential in the squid giant axon. Rubin and Wechselberger [1] reduced this model to a threedimensional vector field in order to understand the role of slow manifolds and canard orbits in organizing the dynamics of the model. In certain situations, the reduced models accurately capture the qualitative behavior of the original higher-dimensional models. However, important aspects of the dynamics may be lost after a reduction [51,52]. Thus, it is very useful to be able to investigate a given system without applying reduction techniques first.
As a natural step for studying higher-dimensional models, we focus on systems with two fast and two slow variables. In particular, we aim to compute and visualize slow manifolds and associated canard orbits in four-dimensional vector fields without the need for any reduction. To this end, we introduce a method for computing twodimensional saddle slow manifolds in R 4 . Computation of saddle slow manifolds is more challenging than that of attracting and repelling slow manifolds due to the existence of both contracting and expanding fast normal directions. This challenge can be overcome with a two-point boundary-value problem (2PBVP) setup. Namely, we compute two-dimensional submanifolds of saddle slow manifolds and their stable and unstable manifolds as families of orbits segments. An important difference with previous studies, is that, in systems with two fast and two slow variables, canard orbits are no longer intersections of slow manifolds. In this paper, we develop and present a homotopy method for computing canard orbits as orbits segments on the attracting slow manifold that stay close to the saddle slow manifolds for a significant amount of time.
We study two different four-dimensional vector fields that feature two-dimensional saddle slow manifolds. The first is an extended version of the normal form of a folded node [34,47], where a fast variable with a stable fast direction is added. In this system, we implement and test our 2PBVP-based approach for computing saddle slow manifolds, and their stable and unstable fast manifolds, as well as associated canard orbits. We show that our methods give reliable approximations of slow manifolds and canard orbits. As a second example, we consider the four-dimensional Hodgkin-Huxley model [8]. Its dynamics is very complex and features periodic MMOs with a series of action potentials separated by subthreshold oscillations. Rubin and Wechselberger [42] nondimensionalized the Hodgkin-Huxley model and obtained a system with two fast and two slow variables, and a simple explicit time-scale separation parameter ε. However, they then focused on a further reduced three-dimensional version of this model. Here, we compute attracting and saddle slow manifolds, as well as associated canard orbits directly in the nondimensionalized four-dimensional Hodgkin-Huxley model. Indeed, we demonstrate that such computations can be done without the need to reduce the system first. Our results confirm and illustrate how slow manifolds and their associated canard orbits play a central role in organizing MMOs also in the four-dimensional model.
The outline of the paper is as follows. Section 2 reviews ideas and concepts from GSPT for systems with two fast and two slow variables. Section 3 begins with a 2PBVP setup for computing orbit segments in four-dimensional vector fields. This is followed by a general method for computing two-dimensional saddle slow manifolds, as well as subsets of their three-dimensional stable and unstable manifolds in systems with two fast and two slow variables. Our approach is then demonstrated and implemented for an extended four-dimensional normal form of a folded node. In Sect. 4, we examine the interaction between attracting and saddle slow manifolds. We also introduce a general approach for computing canard orbits in R 4 and then test it for the extended normal form of a folded node. In Sect. 5, we investigate the underlying dynamics of the full four-dimensional Hodgkin-Huxley model. A bifurcation analysis of MMOs and a short review of GSPT for this model are presented first. Thereafter, we compute attracting and saddle slow manifolds and study their interactions. For a relatively large ε, we illustrate the geometry of so-called ribbons and their bounding twin canard orbits, as well as their role in organizing MMOs. Finally, we perform a continuation analysis of canard orbits of the model. Conclusions, a summary of results and directions for future work can be found in Sect. 6.

Background on Geometric Singular Perturbation Theory
In this paper, we focus on four-dimensional slow-fast (singularly perturbed) systems that can be written in the form with fast variables x ∈ R 2 , slow variables y ∈ R 2 , parameter vector λ ∈ R p , a small parameter 0 < ε 1 representing the time-scale ratio, and sufficiently smooth functions f : The overdot denotes the derivative with respect to the fast time scale t. Geometric singular perturbation theory (GSPT) exploits the separation of different time scales in order to explain the complex dynamics of slow-fast systems. The idea behind GSPT is to analyze two lower-dimensional subsystems of the singular limit and gather information from both subsystems to understand the behavior of the full system. Taking the singular limit of (1) yields the two-dimensional fast subsystem (layer problem) where the slow variables y can be treated as bifurcation parameters. The flow of the fast subsystem is called the fast flow. In contrast, the reduced system (slow subsystem) is obtained by rescaling the time t of (1) to τ = εt and then taking the singular limit. Hence, the reduced system is the two-dimensional system of differential algebraic equations where the prime denotes the derivative with respect to the slow time τ . The flow of the reduced system (3) is called the slow flow. Equilibrium points of the fast subsystem (2) define a two-dimensional manifold called the critical manifold which is also the phase space of the reduced system (3). A submanifold S m ⊆ S is normally hyperbolic if and only if all points on S m are hyperbolic equilibria with respect to the fast variables; that is, the eigenvalues of the Jacobian matrix D x f (x, y, λ, 0) at these points have no zero real part. Depending on these eigenvalues, a corresponding submanifold S m ⊆ S can be attracting (S a ), repelling (S r ) or of saddle type (S s ). According to Fenichel theory, for ε sufficiently small, normally hyperbolic submanifolds of the critical manifold persist as locally invariant slow manifolds [11,12]. A slow manifold S m ε has the same smoothness and stability properties as the corresponding submanifold S m ⊆ S. Moreover, the stable and unstable manifolds, W s (S m ) and W u (S m ), of a normally hyperbolic submanifold S m persist as locally invariant stable and unstable (fast) manifolds, W s (S m ε ) and W u (S m ε ), of the perturbed slow manifold S m ε , respectively. We also remark that the flow on the perturbed slow manifolds is diffeomorphic to the associated slow flow on the critical manifold for sufficiently small ε.
We are interested in the case where an attracting sheet S a meets a saddle sheet S s at a one-dimensional fold curve F . In this situation, the reduced system is singular (not normally hyperbolic) at F . Nevertheless, one attains a description of the dynamics near F by rescaling (3) to obtain the desingularized reduced system where (x, y) ∈ S, and the prime now denotes the derivative with respect to Here, adj(D x f (x, y, λ)) and det(D x f (x, y, λ)) are the adjoint and determinant of the Jacobian matrix D x f (x, y, λ), respectively. Note that this rescaling causes a reversal of the flow on the saddle sheet S s . Equilibrium points of the desingularized reduced system (5) that lie on fold curves are called folded singularities; they can be nodes, saddles or foci. In particular, the existence of a folded node on F allows for a whole sector of trajectories of (3) to pass from S a to S s in the singular limit ε = 0. For ε sufficiently small, this sector gives rise to a finite number of canard orbits, i.e., trajectories that stay close to the attracting slow manifold S a ε and saddle slow manifold S s ε for O(1) time on the slow time scale [32][33][34].

Computing Two-Dimensional Saddle Slow Manifolds
In this section, we introduce a general numerical framework for computing twodimensional saddle slow manifolds and their stable and unstable manifolds in the context of systems with two fast and two slow variables. We also introduce a general method for computing and continuing canard orbits in this setting. We then illustrate our methods throughout with an extended version of the normal form of a folded node.

Two-Point Boundary-Value Problem Setup for Four-Dimensional Slow-Fast Systems
To begin, we introduce a general setup for representing relevant parts of invariant manifolds in R 4 as families of orbit segments by defining suitable 2PBVPs. Throughout the paper, we use this setup extensively to compute and continue invariant objects in systems of the form of (1). The first step for computing orbit segments as solutions of 2PBVPs is to rescale the vector field (1) aṡ where u ∈ R 4 is the state, λ ∈ R p is the parameter vector and H : R 2 × R 2 × R p × R → R 4 is the right-hand side of (1). The parameter T is the total integration time of the computed orbit segment and it is treated as a free parameter; the overdot now represents the derivative with respect to the rescaled time s = t/T . Hence, the orbit segment u(s) is always defined over the interval s ∈ [0, 1]. We define codimension-i and codimension-j hypersurfaces, Ξ 0 and Ξ 1 , where i + j = 4 and then impose the boundary conditions Equation (6) with boundary conditions given by (7) yields a well-posed 2PBVP that defines a one-parameter solution family of orbit segments representing part of a twodimensional manifold of interest. Furthermore, the solution family of this 2PBVP can be computed reliably by employing pseudo-arclength continuation in combination with a 2PBVP solver, such as collocation. A broader overview of this general 2PBVP setup can be found in [48]. Perhaps the simplest case of (7) is when i = 4 and j = 0, such that the corresponding 2PBVP effectively defines an initial value problem. For the case when i = 3 and j = 1, the 2PBVP yields orbit segments that start from a one-dimensional curve and end at a codimension-one submanifold. This case is commonly used to compute twodimensional invariant manifolds as one-parameter families of orbit segments [48]. For systems of dimension higher than three, one-parameter families of orbit segments can also be computed by starting from a two-dimensional surface and ending at another two-dimensional surface. Implementation of the latter case, when i = 2 and j = 2, will be used to compute two-dimensional saddle slow manifolds.

General Approach for Computing Two-Dimensional Slow Manifolds in R 4
Computations of one-dimensional saddle slow manifolds and associated stable and unstable manifolds were performed with a collocation method and numerical integration [24], iterative methods [63] and two-point boundary-value problem (2PBVP) setup [21]. The flow on a one-dimensional saddle slow manifold is very simple because it has only one degree of freedom. In contrast, the flow on a two-dimensional saddle slow manifold may be quite intricate with subsets of trajectories that have different properties. We now present a general approach for computing suitable families of orbit segments to represent the relevant parts of two-dimensional attracting and saddle slow manifolds in R 4 . We consider a system of the form (1). We further assume that the critical manifold S of (1) has a two-dimensional attracting sheet S a and a two-dimensional saddle sheet S s that meet at a generic fold curve F . As was described in Sect. 3.1, we first rescale the vector field (1) to obtain (6), where the parameter T is the integration time associated with (1). We then implement a 2PBVP with the boundary conditions given in (7), where Ξ 0 and Ξ 1 are codimension-i and codimension-j hypersurfaces with i + j = 4. The choice of Ξ 0 and Ξ 1 determines a one-parameter family of solutions that represent a given two-dimensional invariant manifold.

Computing Attracting and Repelling Slow Manifolds in R 4
We compute the attracting slow manifold S a ε in systems of the form (1) by using the boundary conditions where L a is a one-dimensional curve on S a sufficiently away from F , and Σ a is a codimension-one submanifold transverse to the flow. The computation of a repelling slow manifold can be performed in the same way by reversing time. This approach of computing attracting and repelling slow manifolds is a direct generalization to higher dimensions of the method used in [35,38,45,48]; it can also be used for computing attracting and repelling slow manifolds in any system with two slow and an arbitrary number n ≥ 1 of fast variables.

Computing Saddle Slow Manifolds in R 4
Now we introduce a method for computing two-dimensional saddle slow manifolds in R 4 . This is a challenge because one needs to take into account the nature of the flow on the manifold. In addition, the existence of expanding and contracting fast normal directions means that a saddle slow manifold is not uniformly attracting in either forward or backward time. To deal with these issues, we define suitable 2PBVPs to compute different parts of a given saddle slow manifold. To set up the method, one needs some information about the behavior of the slow flow of the reduced system, as well as the stable and unstable eigendirections associated with the saddle sheet of the critical manifold in the singular limit. We aim to represent the saddle slow manifold S s ε in a system of the form (1) as a collection of one-parameter families of orbit segments that approach S s ε along its stable manifold W s (S s ε ) and leave via its unstable manifold W u (S s ε ). This can be achieved by implementing suitable 2PBVPs. Specifically, we represent each part of S s ε as a family of orbit segments satisfying the boundary conditions: Here, Σ 0 and Σ 1 are suitably chosen two-dimensional (codimension-two) sections. In our approach, each section Σ i is defined as the intersection of two threedimensional (codimension-one) submanifolds: Here, Σ 0 and Σ 1 are three-dimensional submanifolds that are chosen transverse to the stable manifold W s (S s ε ) and unstable manifold W u (S s ε ) of the saddle slow manifold S s ε , respectively; their definitions are informed by the (local) eigendirections of the equilibria of the fast subsystem that lie on S s . In contrast, the three-dimensional submanifolds Σ 0 and Σ 1 are constraints that are tailored to the desired submanifold of S s ε to be computed. In this way, we can represent a suitable submanifold of S s ε as a one-parameter family of orbit segments that enter along W s (S s ε ), and then stay extremely close to S s ε before exiting along W u (S s ε ).

Extended Normal Form of a Folded Node
Our approach of computing two-dimensional saddle slow manifolds is general, but requires some knowledge of the slow flow on the critical manifold. To illustrate it, we now introduce a concrete example of a four-dimensional vector field to provide a case study. The normal form of a folded node is a system with one fast and two slow variables [34,47]. Here, we introduce and consider the extended vector field given by which includes a second fast variable w representing a stable fast direction. Since w does not play a role in the dynamics of the other three variables, the dynamics of x, y and z of the original model, including the slow flow on the critical manifold, remains unchanged. Hence, this system is a good test-case model for computing twodimensional saddle slow manifolds in R 4 . The critical manifold of (11) is given by the two-dimensional parabolic surface S := (x, y, z, w) | x = −z 2 , and z = w , (12) which has an attracting sheet S a and a saddle sheet S s that meet at a nondegenerate fold curve F . More precisely, S is composed of The reduced system of (11) is given by  (11) with the phase portrait of the reduced system (14) with μ = 9.2, projected onto the (x, y, z)-space. The lower sheet S a of the critical manifold (gray) is attracting and the upper sheet S s is of saddle type. The fold curve F (gray line) separates the two sheets. The black dot at the origin is the folded node, and the magenta curve ξ s is the strong singular canard It defines the slow flow, where the prime denotes the derivative with respect to the slow time scale τ = εt. The slow flow on the saddle sheet S s is the image of the slow flow on S a under the symmetry The origin is a folded singularity, which is an equilibrium of the desingularized reduced system on F but not of the full system. Its eigenvalues as an equilibrium of the corresponding desingularized reduced system are −μ and −1; thus, the folded singularity is a node for μ > 0. In the calculations that follow, we first consider μ = 9.2. Figure 1 illustrates the slow flow in the neighborhood of the folded node (black dot) on the critical manifold (gray surface), which is projected onto the (x, y, z)space. The fold curve F (gray line) separates the attracting sheet S a and the saddle sheet S s of the critical manifold. The magenta trajectory is the strong singular canard ξ s , which divides S a into a funnel region and a jump region; see [34] for details. There are three types of trajectories on S a and S s that behave differently. Orange and green trajectories on S a , which lie in the funnel region between ξ s and F , move to S s through the folded node in finite time. The black trajectories on S a move towards F and then jump in the fast direction. The black and green trajectories on S s move away from F all the way to infinity. On the other hand, orange trajectories on S s leave from the vicinity of the folded node and make a return to the fold curve F , where they jump in the fast direction.

Slow Manifolds of the Extended Normal Form
Here, we illustrate our approach of computing attracting and saddle slow manifolds with system (11). Since the flow on the attracting slow manifold S a ε and the saddle slow manifold S s ε , for ε small enough, is diffeomorphic to that on S a and S s , respectively, we use the nature of the slow flow on the respective parts of the critical manifold shown in Fig. 1 as a guide for selecting families of orbit segments that represent the slow manifolds. With this selection, we compute the slow manifolds with the general approach introduced in Sect. 3.2.

Computing the Saddle Slow Manifold of the Normal Form
We represent the saddle slow manifold S s ε of (11) as families of orbit segments that satisfy the boundary conditions given by (9) and (10). The three-dimensional submanifolds Σ 0 and Σ 1 of (10) need to be transverse to the stable manifold W s (S s ε ) and unstable manifold W u (S s ε ) of the saddle slow manifold S s ε , respectively. The extended normal form has the special property that each equilibrium of the fast subsystem on S s has the w-direction as its stable eigendirection, whereas the unstable eigendirection is defined by (x, z) = (2z + 1, 1). This property implies that, away from F , trajectories on W s (S s ε ) always enter S s ε via the w-nullcline, and trajectories on W u (S s ε ) always move away from S s ε via the z-nullcline. Hence, we choose which guarantees that Σ 0 and Σ 1 are transverse to W s (S s ε ) and W u (S s ε ), respectively. We use these choices for Σ 0 and Σ 1 for all computed submanifolds of S s ε . On the other hand, sections Σ 0 and Σ 1 determine the constraints that specify the computed submanifolds of S s ε . They are different for each of the three families of orbit segments shown in Fig. 1. More specifically, the submanifold of S s ε corresponding to the black orbit segments in Fig. 1 is specified by so that we capture the start of these orbit segments near the fold curve F and follow them up to z = 30. The submanifold of S s ε corresponding to the green orbit segments is specified by because these orbit segments all start at y = 0 near the folded node. Finally, the submanifold of S s ε corresponding to the orange orbit segments is specified by because these orbits segments end close to F . The 2PBVPs of the three submanifolds that constitute the computed saddle slow manifold are defined by Eq. (6) with boundary conditions (9), (10) as specified by (16), and one of (17), (18) or (19). Figure 2 shows the three computed submanifolds of the saddle slow manifold S s ε for ε = 0.01. They correspond to the three different families of trajectories on S s shown in Fig. 1. The left column shows each submanifold, projected onto the (x, y, z)-space. The sections Σ 0 (light blue) and Σ 1 (dark blue) used for the boundary conditions are also shown. The right column of Fig. 2 shows the same submanifolds of S s ε , projected onto the (x, z)-and (w, z)-planes. In these projections, the Three computed submanifolds of S s ε of system (11) with μ = 9.2 and ε = 0.01. The left column illustrates the three submanifolds, projected onto the (x, y, z)-space, together with the corresponding boundary conditions Σ 0 and Σ 1 (light-and dark-blue sections). The right column shows the three submanifolds, projected onto the (x, z)-plane and the (w, z)-plane, where we also include stable (blue) and unstable (red) eigenvector directions at a selection of points on S s saddle branch S s of the critical manifold correspond to {(x, z) | x = −z 2 and z > 0} and {(w, z) | z = w and z > 0}, respectively. In both projections, it is notable that all computed orbit segments of S s ε stay extremely close to S s ; this is an indication that they give a good approximation of the saddle slow manifold. In the (w, z)-plane, the short blue (red) segments are the directions of the stable (unstable) eigenvectors associated with the equilibrium points of S s at the singular limit ε = 0. Figure 2(a) shows the computed orbit segments (black) of the first family of S s ε . The green surface, which is a part of S s ε , has been rendered from these black orbit segments. Indeed, the black family of orbit segments are topologically as the black trajectories on S s shown in Fig. 1. In Fig. 2(a1), the sections Σ 0 (light-blue plane) and Σ 1 (dark-blue line) correspond to the boundary conditions that define the computed black orbit segments. The section Σ 1 is actually two dimensional but it appears as a one-dimensional line in projection onto the (x, y, z)-space because the corresponding equation does not depend on w. The boundary conditions for the shown submanifold of S s ε are defined by (9), (10), (16), and (17). As a result, the computed family of orbit segments start from the vicinity of the fold curve F and end far away from F . Panels (a2) and (a3) of Fig. 2 show that this family of S s ε is very close to S s . Observe from Fig. 2(a3) that the stable eigenvectors are aligned with the w-axis. Also note that, sufficiently far away from F , the unstable eigenvectors are (almost) aligned with the z-axis. These eigendirections illustrate why we force the computed orbit segments representing this part of S s ε to start from the w-nullcline and end at the z-nullcline; see (16). In this way, we compute orbit segments of S s ε without including the fast segments of the associated stable and unstable manifolds. However, as we will show in Sect. 3.5, we are able to change the boundary conditions and extend the orbit segments, in backward and forward time, to include fast segments of W s (S s ε ) and W u (S s ε ), respectively. Figure 2(b) shows the computed orbit segments (green) of the second family of S s ε . This family is topologically as the green family of trajectories on S s shown in Fig. 1. The green surface, which is a part of S s ε , has been rendered from the shown green orbit segments. Sections Σ 0 (light-blue plane) and Σ 1 (dark-blue line) illustrate the respective boundary conditions. Here, the section Σ 1 is also two dimensional but it again appears as a one-dimensional curve in projection onto the (x, y, z)-space. The 2PBVP setup for this submanifold of S s ε is defined by (9), (10), (16), and (18). The green orbit segments start from the neighborhood of the folded node and terminate sufficiently far away from the fold curve. In both projections shown in panels (b2) and (b3), the computed orbit segments lie extremely close to S s , which indicates that the computed orbit segments again provide a good approximation of the saddle slow manifold.
Finally, Fig. 2(c) shows the computed orbit segments (orange) of the third family of S s ε . The boundary conditions for this submanifold are (9), (10), (16) and (19). The corresponding orbits segments start from the neighborhood of the folded node and make a return to the fold curve F . Panel (c3) shows that the orbit segments divert slightly from the critical manifold S as they return to the vicinity of F . This is because the linearization of the unstable manifold W u (S s ), which is given by (z, w) = (1 + 2z, 1), becomes more aligned with the slow flow on the critical manifold as F is approached. In other words, the normal hyperbolicity of the saddle slow manifold is lost near F . Nevertheless, the projections in panels (c2) and (c3) indicate that the computed orbit segments still provide a good approximation of the saddle slow manifold.

The Overall Geometry of Slow Manifolds
As was described in Sect. 3.2, we compute the attracting slow manifold S a ε by starting from a line away from the fold curve F and ending at a three-dimensional submanifold. More specifically, we represent S a ε as families of orbit segments that are subject to the boundary condition given by (8), where the precise loci of L a and Σ a depend on the desired submanifold of S a ε to be computed. To represent the black, green and orange trajectories of S a in Fig. 1 as orbit segments of S a ε , we choose and respectively. These boundary conditions result in three families of orbit segments that correspond to three submanifolds constituting the attracting slow manifold S a ε . Figure 3 shows the resulting approximations of the attracting slow manifold S a ε (red surface) and the saddle slow manifold S s ε (green surface) for ε = 0.01; also shown are the computed orbit segments used to render the surfaces. Note that the flow on the slow manifolds is qualitatively the same as the slow flow of the reduced system (ε = 0) shown in Fig. 1, as one may expect for such a small value of ε = 0.01. The saddle slow manifold S s ε is rendered as a concatenation of the three families of orbit segments as shown in Fig. 2. Likewise, the attracting slow manifold S a ε is rendered as a concatenation of the three computed families of orbit segments based on the boundary conditions (20), (21) and (22), respectively. The flow on S a ε is topologically equivalent to the slow flow on S a in Fig. 1. Black trajectories of S a ε approach Fig. 3 Attracting slow manifold S a ε (red surface) and saddle slow manifold S s ε (green surface), projected onto (x, y, z)-space, for system (11) with μ = 9.2 and ε = 0.01. The gray line is the fold curve F of the critical manifold. These surfaces were rendered from the orbit segments shown; compare with Fig. 2 the fold curve F , near which they are expected to make a jump away from the critical manifold. In contrast, green and orange trajectories of S a ε terminate near the folded node where they interact with the saddle slow manifold. The interaction between S a ε and S s ε will be discussed in Sect. 4.

Computing Stable and Unstable Manifolds of Saddle Slow Manifolds
Two-dimensional saddle slow manifolds in R 4 are associated with three-dimensional stable and unstable (fast) manifolds. It is very challenging to compute and visualize such three-dimensional geometric objects. Nevertheless, we now present a general approach for computing two-dimensional submanifolds of the three-dimensional stable manifold W s (S s ε ) and the three-dimensional unstable manifold W u (S s ε ) of the saddle slow manifold S s ε . This can be achieved by implementing suitable 2PBVPs to compute orbit segments that start transversally to W s (S s ε ) and end transversally to W u (S s ε ). The basic idea of our approach is inspired by the computation of twodimensional (un)stable manifolds of a one-dimensional saddle slow manifold, which is presented in [21].
Consider system (1) and assume that it features a two-dimensional saddle slow manifold S s ε with associated stable manifold W s (S s ε ) and unstable manifold W u (S s ε ). Each normally hyperbolic trajectory c s ε ⊂ S s ε on the saddle slow manifold is associated with two-dimensional stable and unstable manifolds, W s (c s ε ) and W u (c s ε ), that approach the (slow) trajectory c s ε at an exponential rate in forward and backward time, respectively. To compute the stable manifold W s (c s ε ) of a trajectory c s ε , we first extend c s ε by the flow in backward time to include a fast stable segment. This can be done for each of the two sides of the corresponding stable eigenvectors. We then continue the resulting extended orbit segment as a one-parameter family of orbit segments subject to the boundary conditions where Σ s is a codimension-one submanifold far away from S and transverse to the stable manifold of the trajectory c s ε . On the other hand, L s is a line segment close to the saddle sheet S s and transverse to the unstable manifold of the trajectory c s ε . We can also compute two-dimensional submanifolds of the stable manifold W s (S s ε ) associated with a one-dimensional curve on S s ε . This can be achieved by imposing the boundary conditions where Σ s 1 and Σ s 0 are codimension-two submanifolds that are suitably chosen to be transverse to W s (S s ε ) and W u (S s ε ), respectively. Submanifolds of the unstable manifold W u (S s ε ) can be computed with the same setup after reversing time. The precise definition of the 2PBVPs given by (23) and (24) depends on the desired computed submanifold, as will be described for our specific example in the next section. Figure 4 shows the stable manifold W s (c s ε ) (blue surfaces) and unstable manifold W u (c s ε ) (red surfaces) of a selected trajectory c s ε (green curve) on the saddle slow manifold S s ε , projected onto the (w, y, z)-space. Figure 4(a) shows an approximation of the two-dimensional local manifold W s (c s ε ) (blue surface); also shown is a selection of orbit segments on W s (c s ε ) (blue curves). Each side of the local manifold is computed by starting from a chosen orbit segment, which is then continued as a one-parameter family that satisfies the boundary conditions Here, Σ s is a codimension-one submanifold and L s is a one-dimensional curve. The submanifold Σ s is chosen to be transverse to W s (c s ε ) and sufficiently far away from S. The curve L s is chosen to be transverse to W u (c s ε ). This computational setup forces the solution segments to move along W s (c s ε ) before converging to the slow trajectory c s ε (in the direction of increasing z). Figure 4(b) shows the local manifold W u (S s ε ) (red surfaces) of the same trajectory c s ε (green curve), together with a selection of orbit segments on W u (c s ε ) (red curves). Again, we compute each side of the local manifold by starting from a chosen orbit segment, which is then continued as a one-parameter family that satisfies the boundary conditions

Stable and Unstable Manifolds of a Trajectory on S s ε
where L u and Σ u are chosen to be transverse to W s (c s ε ) and W u (c s ε ), respectively. The computed orbit segments move along the slow trajectory c s ε (again in the direction of increasing z), before exiting via its two-dimensional unstable manifold W u (c s ε ).

Slices of Stable and Unstable Manifolds of S s ε
To compute different parts of the three-dimensional manifolds W s (S s ε ) and W u (S s ε ), we implement variations of suitable 2PBVPs. Figure 5 shows a two-dimensional piece of S s ε (green surface), together with six two-dimensional submanifolds of the associated three-dimensional stable manifold W s (S s ε ) (blue surfaces) and unstable manifold W u (S s ε ) (red surfaces) of system (11). All surfaces are rendered from the respective computed orbit segments (thick curves). Figure 5(a) shows six two-dimensional submanifolds of W s (S s ε ), each formed by a collection of particular orbit segments on W s (c s ε ) for a family of slow trajectories c s ε on S s ε . More precisely, each submanifold corresponds to a one-parameter family of orbit segments that satisfy the boundary conditions where the sections Σ s 0 and Σ s 1 are two-dimensional constraints that are chosen to be transverse to W s (S s ε ) and W u (S s ε ), respectively, and w s is a suitably chosen constant that identifies a particular selection of orbit segments on W s (S s ε ). The six selected submanifolds of W s (S s ε ) are defined by w s = −19, −10, 0, 20, 30, and 40. Similarily, Fig. 5(b) shows six two-dimensional submanifolds of the threedimensional unstable manifold W u (S s ε ). Each submanifold corresponds to a one-

Fig. 5
A selection of two-dimensional submanifolds of the three-dimensional manifolds W s (S s ε ) and W u (S s ε ) of system (11) with μ = 9.2 and ε = 0.01. Panel (a) shows a piece of S s ε (green surface) and associated two-dimensional submanifolds of W s (S s ε ) (blue surface). Panel (b) shows a piece of S s ε (green surface) and associated two-dimensional submanifolds of W u (S s ε ) (red surface). All surfaces are rendered from the respective thick orbit segments parameter family of orbit segments that satisfy the boundary conditions

Interaction Between Attracting and Saddle Slow Manifolds and Canard Orbits in R 4
Slow manifolds can be extended by the flow up to a transverse section near a folded node where they interact with each other. In systems of ordinary differential equations with one fast and two slow variables, extended attracting and repelling slow manifolds may intersect transversally in isolated canard orbits. In such systems, canard orbits can be detected by placing a cross-section transverse to both manifolds near the folded node and determining their intersections; see also [35,45,48,64]. In systems with two fast and two slow variables, attracting and saddle slow manifolds spiral around each other in forward and backward time, respectively, in the vicinity of the folded node. However, their possible intersections do not occur in a structurally stable manner in R 4 . We now investigate the interaction between two-dimensional attracting and saddle slow manifolds of the four-dimensional system (11). In the computations that follow, we set μ = 100.1 in order to obtain more interesting dynamics near the folded node, that is, more canard orbits near ε = 0. Figure 6 illustrates how the attracting slow manifold S a ε and saddle slow manifold S s ε interact with each other near the folded node (the origin) in system (11) with μ = 100.1 and ε = 0.01. Panel (a) shows a submanifold of S a ε (red surface) and a submanifold of S s ε (green surface), computed up to the three-dimensional section Σ := {y = 0} (blue plane), which contains the folded node at the origin; the figure shows a projection onto the (x, y, z)-space. Each surface intersects Σ in a onedimensional curve. Panel (b) shows the intersection curves S a ε = S a ε ∪ Σ (red) and S s ε = S s ε ∪ Σ (green) in the (x, w, z)-space representing Σ . The two curves spiral around each other but do not actually intersect. Panel (c1) illustrates the same two curves in simultaneous projections onto the (x, z)-and (x, w)-planes; and panel (c2) is an enlargement. The curves S a ε,z (red) and S a ε,w (orange) are the projections of S a ε onto the (x, z)-and (x, w)-planes, respectively. The curve S s ε,z = S s ε,w (green) is the projection of S s ε , which is identical in both projections because it lies on Σ 1 so that z = w. An intersection between the two surfaces occurs when the three projected curves S a ε,z , S a ε,w and S s ε,z = S s ε,w all coincide in a single point. Hence, panels (c1) and (c2) show that there is no transversal intersection of S a ε and S s ε , which represents the generic situation. On the other hand, there are many near intersections, which implies that there should be a number of orbit segments that lie on S a ε , stay extremely close to S s ε for a finite time and then leave via W u (S s ε ). Such orbit segments are canard orbits as we discuss next.

General Approach for Computing Canard Orbits in R 4
In systems with one fast and two slow variables, canard orbits are defined as transversal intersections of slow manifolds [33,34]. In systems with two fast and two slow variables, on the other hand, attracting and saddle slow manifolds are two dimensional and generically do not intersect, as was demonstrated in the previous section. This does not mean, however, that canard orbits do not exist in this higherdimensional setting. Rather, canard orbits are trajectories that follow attracting and saddle slow manifolds for O(1) time on the slow time scale before exiting via the (red) and S s ε (green) projected onto the (x, y, z)-space, and computed up to the three-dimensional section Σ (blue); the gray surface is the critical manifold S and the gray line is the fold curve F . Panel (b) shows the intersection curves S a ε = S a ε ∪ Σ (red) and S s ε = S s ε ∪ Σ (green) in the (x, w, z)-space of Σ . Panel (c1) shows the curves S a ε,z (red) and S a ε,w (orange) as projections of S a ε onto the (x, z)-and (x, w)-planes, respectively, together with the intersection curve S s ε (green), which is the same in both projections. Panel (c2) is an enlargement of panel (c1) unstable manifold of the saddle slow manifold [65]. This well-known alternative definition of a canard orbit has the advantage that it does not depend on the dimension of the slow manifold. As such, we are able to propose a general approach for computing and continuing canard orbits as orbit segments that follow S a ε and stay close to S s ε for a significant amount of time before leaving via W u (S s ε ). Consider system (1) and assume that the critical manifold has attracting and saddle sheets, S a and S s , respectively, that meet at a generic fold curve F on which there is folded node. We compute canard orbits as orbit segments that satisfy the boundary conditions Here, L a is a one-dimensional curve on S a sufficiently far away from F . The section Σ 1 is the two-dimensional intersection of three-dimensional submanifolds Σ 1 and Σ 1 , where Σ 1 is chosen to be transverse to W u (S s ε ), and Σ 1 must be sufficiently far away from F . Note that Σ 1 is chosen in the same way as for computing the saddle slow manifold; see (10). This computational setup forces the computed solution segments to follow S a ε and then S s ε for a sufficient amount of time before leaving via the unstable manifold W u (S s ε ). Such orbit segments provide excellent approximations of canard orbits, in the same spirit as those computed in [38].
To obtain a solution segment that satisfies (29), one needs to perform a homotopy step. The underlying idea is to start from an orbit segment with u(0) ∈ L a and u(1) on either Σ 1 or Σ 1 . We then continue this orbit segment until (29) is satisfied. Once detected, canard orbits that are solutions of the well-posed 2PBVP with the boundary conditions given by (29) can be continued in any system parameter.

Detecting and Continuing Canards of the Normal Form
To detect canard orbits of system (11), we represent them as orbit segments that satisfy (29) with the system-specific choice Note that L a 1 is the same line as was used for computing the first and second submanifolds of S a ε , and Σ 1 is the same section used for computing the first and second submanifolds of S s ε . Nonetheless, solutions of (30) do not correspond to transversal intersections of S a ε and S s ε . To detect a canard orbit ξ i with i small-amplitude oscillations (SAOs), we start from an orbit segment on S a ε that exhibits i SAOs. We then extend this orbit segment by the flow until the end point reaches the section Σ 1 = {z = 30}. We then continue this orbit segment by imposing the boundary conditions and letting u(0) vary along L a 1 . This yields a one-parameter family of orbit segments on S a ε with i SAOs. As soon as u(1) lies on Σ 1 = {ż = 0}, Eq. (30) is satisfied and the corresponding orbit segment represents a canard orbit ξ i . Figure 7 illustrates this homotopy step for finding canard orbits in system (11) with μ = 100.1 and ε = 0.01. Panels (a) and (b) show orbit segments with four SAOs that satisfy conditions (31), in projection onto the (y, z)-plane and (x, z)-plane, respectively. These orbits segments stay very close to S s ε for a certain amount of time before leaving via W u (S s ε ). However, only the red orbit segment satisfies (30) and stays close to S s ε for a longer time than all of the other computed orbit segments. This is because this orbit segment 'switches' to leave via W u (S s ε ) from the other side of the stable manifold W s (S s ε ). Hence, we regard it as the representative canard orbit ξ 4 with four SAOs. This homotopy step is an extension of the approach of detecting canard orbits in R 3 that was presented in [38]. For the parameter value μ = 100.1 and sufficiently small ε, the theory predicts the existence of two primary canards and 50 additional secondary canards [34]. Each of these canard orbits can be detected as solutions that satisfy (30), via the homotopy step discussed in Sect. 4.1.2. Figure 8 shows nine selected canard orbits ξ 0 -ξ 8 of system (11) with μ = 100.1 and ε = 0.01. Panel (a) shows the canard orbits in projection onto the (y, z)-plane, while the top left and bottom right insets show projections onto the (x, z)-and (w, z)-planes, respectively; panel (b) is an enlargement. The primary canard orbit ξ 0 follows both S a ε and S s ε without making any rotation. This canard orbit is the equivalent of the strong singular canard ξ s at the singular limit ε = 0. Each canard orbit ξ i , i > 0, makes i rotations around the singular weak canard ξ w (not shown). Observe that the projections of canard orbits in the insets are close to those of the critical manifold S; see Eq. 4. This provides a strong evidence that all of the computed canard orbits ξ 0 -ξ 8 stay extremely close to both S a and S s . This is a clear indication that the computed orbit segments are accurate approximations of canard orbits.

Continuation of Canard Orbits
In this section, we perform an analysis of how the canard orbits depend on μ and ε. This is readily achieved by continuing them as orbit segments that satisfy the boundary conditions given by (30). Figure 9 shows the continuation of canard orbits ξ 0 -ξ 8 in ε. We again use ξ i to denote the ε-dependent branches that correspond to canard orbits with i SAOs. Note that all branches converge to the singular limit ε = 0. Indeed, we find that all canard orbits accumulate on the strong singular canard ξ s near the singular limit. This agrees with predictions by the theory [34]. Figure 10 shows the continuation of canard orbits ξ 0 -ξ 8 [34], the branches ξ i are expected to terminate exactly at μ = 2i + 1 near the singular limit ε = 0. Again, our numerical results are consistent with the theory, which is another indication of the reliability of our computations.

Slow Manifolds and Canard Orbits in the Full Hodgkin-Huxley Model
In this section, we implement and demonstrate our numerical approach for computing two-dimensional saddle slow manifolds and associated canard orbits in an exam- Fig. 9 Continuation in ε of canard orbits ξ 0 -ξ 8 for system (11) with μ = 100.1

Fig. 10
Continuation in μ of canard orbits ξ 0 -ξ 8 for system (11) with ε = 0.01 ple from neuroscience, namely, the Hodgkin-Huxley model. The Hodgkin-Huxley model was originally formulated to describe the action potential of the squid giant axon [8]. We study a singularly perturbed nondimensionalized version of the Hodgkin-Huxley model obtained by Rubin and Wechselberger [1], given by where v represents the action potential and m, h and n represent gating variables. The parameter I > 0 is the injected current and ε is the time-scale separation parameter. The steady-state functions m ∞ (v), h ∞ (v) and n ∞ (v), and time functions t m (v), t h (v) and t m (v) are standard from the Hodgkin-Huxley formalism, and they are given in the Appendix. In this paper, we vary I and ε and keep the other parameters fixed as in Table 1; the given parameter values are the same used in [1] and they correspond to the original Hodgkin-Huxley model, except that we set τ h = 2 instead of τ h = 1 to obtain more interesting dynamics. We initially set I = 9.74 and ε = 0.0083. For ε small, v and m can be regarded as fast variables, and h and n as slow variables. Rubin and Wechselberger [1] applied a center manifold reduction to (32) by setting m = m ∞ (v) to eliminate m and obtain a three-dimensional reduced model. They applied methods from GSPT to prove the existence of relaxation oscillations, MMOs and canard orbits in the three-dimensional reduced model. They also performed bifurcation analyses in I and ε for various choices of τ h and τ n [42]. Desroches et al. [36] detected and continued canard orbits for the same reduced three-dimensional model. Other versions of the full four-dimensional Hodgkin-Huxley model have been investigated in previous studies [8,[66][67][68][69]. Here, we study the nondimensionalized version of the full Hodgkin-Huxley model (32) without applying any reduction. More specifically, we examine the role of slow manifolds and associated canard orbits in organizing MMOs.

Bifurcation Diagram of the Four-Dimensional Hodgkin-Huxley Model
We start our analysis with the bifurcation diagram of the Hodgkin-Huxley model (32) in the parameter I for fixed ε = 0.0083. Figure 11 shows the L 2 -norm of equilibria (black curve) and periodic orbits (colored curves) versus I ; solid and dashed curves indicate stable and unstable branches, respectively. Note that there is a single equilibrium point for I ∈ [0, 200]. It undergoes subcritical and supercritical Hopf bifurcations (black dots) at I ≈ 9.70009 and I ≈ 170.609 (labeled HB 1 and HB 2 ), respectively. In panel (a), the primary branch of periodic orbits (blue curve) connects HB 1 with HB 2 . Panel (b) shows an enlargement for I ∈ [8.5986, 14.5045]; stability properties are not shown here. The primary branch undergoes a period-doubling bifurcation (PD) at I ≈ 8.80514. The period-doubled branch (cyan) that emanates from PD terminates at a branch point (BP) at I ≈ 14.4103. Figure 11(b) also shows 17 other branches of periodic orbits that are alternately colored gray, red and green; these branches are isolas and some of them have folds extremely close to the point BP. Panel (c) shows an enlargement of panel (b) with stability properties shown. Each branch goes through a number of saddle-node bifurcations of periodic orbits and period-doubling bifurcations. Therefore, different types of stabilities can be found along these branches. The labels indicate the MMO signatures of the stable parts at  [19,70]. The periodic orbits that lie on a given branch are topologically equivalent; however, the oscillations exhibited by the periodic orbits vary in amplitude along that branch. For example, the branch labeled 5 1 also includes periodic orbits with signatures of type L 1 It is notable that some branches coexist for a range of parameter values. For instance, we find at least 26 coexisting periodic orbits for I = 9.74. For the remainder of this paper, we fix I = 9.74 and investigate the structure of the slow manifolds and other dynamics of (32).

Slow-Fast Analysis of the Four-Dimensional Hodgkin-Huxley Model
To gain insight into the geometric mechanism of MMOs, we start with applying methods from GSPT to the four-dimensional Hodgkin-Huxley model. Taking the singular limit ε = 0 of (32) yields the two-dimensional fast subsystem v, m, h, n), where h and n can be treated as bifurcation parameters. Equilibrium points of (33) form a two-dimensional surface corresponding to the critical manifold S, which is a cubic-shaped smooth surface with two fold curves. Figure 12 shows the critical manifold for I = 9.74 in two different projections. The two fold curves F 1 and F 2 are generic fold (saddle-node) bifurcations of (33), and they separate the critical manifold into three different sheets. Note that the values of v and m hardly change along the upper fold curve F 2 . The lower and upper sheets (S a , S a ) of the critical manifold are attracting since the corresponding two normal eigenvalues of the equilibria are negative. The middle sheet S s is of saddle type since the corresponding equilibria have eigenvalues of opposite signs. The reduced system of (32) is the two-dimensional system of differential algebraic equations: m, h, n), where the prime denotes the derivative with respect to the slow time scale τ = εt.
With the implicit function theorem, the reduced system can be written in the form Straightforward calculations show that system (35) is singular along the fold curves F 1 and F 2 . This problem can be avoided via time rescaling by the factor −( ∂f where the prime now denotes the derivative with respect to the new rescaled time System (36) is the desingularized reduced system, for which the orientation on the saddle sheet S s of the critical manifold is reversed. We found that, for I = 9.74, the desingularized reduced system (36) features two equilibrium points. One equilibrium lies on the saddle sheet S s and corresponds to a saddle-focus equilibrium point q of the full system (32), given by q = (v, m, h, n) ≈ (−0.596423, 0.973836, 0.405820, 0.401974).
This equilibrium point has a two-dimensional stable manifold W s (q) and a twodimensional unstable manifold W u (q). The other equilibrium of (36) lies on the lower fold curve F 1 and corresponds to the folded node singularity: The strong and weak stable manifolds of p in system (36) correspond to the strong singular canard ξ s and the weak singular canard ξ w . This allows for the existence of the two-dimensional funnel region, on which solutions of system (35) move from the lower attracting sheet S a through p to the saddle sheet S s in finite time; compare with Fig. 1. The eigenvalues of p are approximately −0.003945 and −0.000125 and their ratio is 31.56. Based on this ratio, the theory [34] predicts that the folded node gives rise to two primary canard orbits and 15 secondary canard orbit near the singular limit ε = 0.
Recall that, for ε sufficiently small, the sheets S a and S s perturb to attracting and saddle slow manifolds S a ε and S s ε , respectively. Additionally, the stable W s (S s ) and unstable manifolds W u (S s ) of the saddle sheet S s also persist as stable W s (S s ε ) and unstable W u (S s ε ) fast manifolds of the saddle slow manifold S s ε . We aim to use the numerical techniques introduced in Sect. 3 to compute slow manifolds and associated canard orbits, and subsequently examine the local and global influence of their geometry on the dynamics of the Hodgkin-Huxley model (32).

Computing the Attracting Slow Manifold
In order to understand the underlying dynamics of the Hodgkin-Huxley model (32), we perform computations of the locally invariant slow manifolds. The attracting slow manifold S a ε of system (32) can be computed with the generalized setup for computing attracting and repelling slow manifolds that was introduced in Sect. 3.2. More specifically, we represent submanifolds of S a ε as families of orbit segments that satisfy the 2PBVP setup given by (8), where L a is a line on the lower attracting sheet S a sufficiently far away from the lower fold curve F 1 , and Σ a is a suitably chosen threedimensional submanifold transverse to the attracting slow manifold. Figure 13 shows two computed submanifolds of the attracting slow manifold S a ε associated with the lower attracting sheet S a . Panel (a) shows the first submanifold of S a ε (red surface), which is rendered from a family of orbit segments (black) that satisfy the systemspecific boundary conditions Figure 13(b) shows the second submanifold of S a ε (red surface), which is rendered with a different family of orbit segments (green), satisfying the boundary conditions Panel (c) shows the overall attracting slow manifold S a ε near p consisting of the two computed families of orbit segments shown in panels (a) and (b).

Computing the Saddle Slow Manifold
We now implement the approach introduced in Sect. 3.2 to compute the twodimensional saddle slow manifold S s ε of (32) by representing it as a family of orbit segments computed with the 2PBVP setup given by (9) and (10). As for S a ε , we represent the saddle slow manifold S s ε with two submanifolds, which are shown in Fig. 14. Panel (a) shows a green surface, which is rendered from the first family of orbit segments (green). Panel (b) shows the second submanifold of S s ε (green surface) rendered from a different family of orbit segments (black). In panels (a) and (b), the sections Σ 0 and Σ 1 (blue planes) correspond to the boundary conditions of the respective computed orbit segments. These sections are defined by (10), where the  Panels (a) and (b) show two submanifolds of S a ε (red surfaces) together with their respective boundary conditions L a 1 and L a 2 (blue lines), and Σ a (blue planes); they are rendered from the orbit segments shown as thick curves (black and green, respectively). Panel (c) shows both submanifolds of S a ε together as a single surface three-dimensional submanifolds Σ 0 , Σ 1 , Σ 0 and Σ 1 are chosen as follows. The sections Σ 0 and Σ 1 of (10) must be transverse to the associated stable manifold W s (S s ε ) and unstable manifold W u (S s ε ), respectively. We use a different section Σ 0 for each submanifold of the computed saddle slow manifold to ensure that the corresponding orbit segments are transverse to W s (S s ε ). Specifically, we choose for the black orbit segments, and  Panels (a) and (b) show two submanifolds of S s ε (green surfaces) together with their respective boundary conditions Σ 0 and Σ 1 (blue planes); they are rendered from the orbit segments shown as thick curves (black and green, respectively). Panel (c) shows both submanifolds of S s ε together as a single surface for the green orbit segments. The section is used for both submanifolds of S s ε , and it is chosen transverse to W u (S s ε ). In contrast, sections Σ 0 and Σ 1 determine the constraints of the computed submanifold of S s ε . More specifically, for the green orbit segments, we choose so that these orbit segments start near the fold curve F 1 and end at v = −0.28. For the black orbit segments, we choose which means that these orbit segments start from h = 0.55 instead. Figure 14(c) shows both submanifolds together, so that they form the entire surface, that is, the saddle slow manifold S s ε .

Interaction Between Attracting and Saddle Slow Manifolds
Attracting and saddle slow manifolds of system (32) can be extended by the flow in forward and backward time, respectively, to study the SAOs that occur in the vicinity of the lower fold curve F 1 . The mechanism for these SAOs can be explained by the combined presence of the folded node p given by (39) and the saddle-focus equilibrium q given by (38). Figure 15 shows the computations of the attracting slow manifold S a ε and the saddle slow manifold S s ε extended up to the three-dimensional cross section Σ := {h = 0.4}. Panel (a) shows S a ε (red) and S s ε (green) in projection onto (n, h, v)-space. Panel (b) shows the intersection curves S a ε ∩ Σ (red) and S s ε ∩ Σ (green) in the (m, n, v)-space representing Σ . Since they are again very close to each other, it is very challenging to know whether there are any intersections between the two curves. Therefore, we use a different projection technique to gain a better insight. Panel (c) shows the intersection curves S a ε ∩ Σ and S s ε ∩ Σ , projected onto the (n, v)-plane. Panel (d) illustrates the same intersection curves with the values of m color coded as shown by the color bar; panels (e) and (f) are enlargements of panel (d). We find no intersections between the curves S a ε ∩ Σ and S s ε ∩ Σ , as is expected for two curves in R 4 . Nevertheless, the colors show that there are again many near intersections. Hence, as before, there should be orbit segments on S a ε that stay extremely close to the saddle slow manifold S s ε before leaving via its unstable manifold W u (S s ε ). Such orbit segments are canard orbits in R 4 , as was discussed Sect. 4.

Computing Canard Orbits of the Hodgkin-Huxley Model
We now implement the general approach introduced in Sect. 4.1 to detect canard orbits of system (32). More specifically, we compute them as orbit segments that lie on the attracting slow manifold S a ε and stay extremely close to the saddle slow manifold S s ε before leaving via its unstable manifold W u (S s ε ); see equations (29). Here, we use a different homotopy step from the one used for system (11). To detect a canard orbit ξ i with i SAOs, we start from a solution that lies on S a ε and exhibits i SAOs in the vicinity of the lower fold F 1 . We then continue this solution as a one-parameter family of orbit segments subject to the boundary conditions Here, L a 2 is the line defined in (41) that we used for computing the second submanifold of S a ε , and Σ 1 is the section defined by (44). These boundary conditions (47) yield a family of orbit segments on S a ε that terminate at the v-nullcline, which is transverse to the unstable manifold W u (S s ε ). They also guarantee that the orbit segments exhibit exactly i SAOs. Figure 16 shows this homotopy step for computing canard orbits for system (32). Panel (a) illustrates how the integration time T depends on the v-coordinate v 1 of the end points on Σ 1 for a family of orbit segments that satisfy (47); the black, gray and red dots correspond to a selection of the computed orbits segments. Panel (b) shows the orbit segments that correspond to the dots of panel (a) in their respective colors, in projection onto the (n, v)-plane. The shown orbit segments start from a very small interval on L a 2 . More specifically, their initial h-values along L a 2 agree to 6 decimal places, namely, h ≈ 0.135841. All shown orbits follow S a ε , exhibit four SAOs, and then follow S s ε for a certain amount of time before leaving via W u (S s ε ). Thus, these orbit segments are indeed transversal intersections of S a ε and W u (S s ε ). The thick red orbit segment corresponds to the trajectory that stays close to S s ε for the longest time, and as such, we regard it as the representative canard orbit ξ 4 .
We performed the same homotopy step to detect other canard orbits ξ i with i SAOs in system (32) and found that, for all i, the orbit segment that spends the longest time near S s ε has a v-maximum of v ≈ −0.28 rounded to the second decimal place. Therefore, we compute all canard orbits as orbit segments that satisfy (47) and have a v-maximum of v = −0.28. More precisely, we compute all canard orbits of system (32) as orbit segments that satisfy the boundary conditions Note that Σ 1 is the section we used for computing both submanifolds of S s ε ; see equations (44), (45) and (46).

Selection of Canard
Orbits for I = 9.74 Figure 17 shows six selected canard orbits (thick curves) together with the extended attracting and saddle slow manifolds. The slow manifolds S a ε (red surface) and S s ε (green surface) are extended by the flow to include the SAOs in the vicinity of the lower fold curve F 1 . The primary canard orbit ξ 0 separates trajectories on S a ε that make at least one rotation around F from those that escape the fold region without making any rotations. Each canard orbit ξ i , i > 0, makes i SAOs in the vicinity of the fold. Recall that there should exist 15 secondary canard orbits for I = 9.74 near the singular limit ε = 0. Nevertheless, we detect more than 111 canard orbits for ε = 0.0083, and ξ 111 is shown as part of the selection in Fig. 17.

Ribbons of the Attracting Slow Manifolds
In [38], we showed that the signature of a given MMO in the three-dimensional autocatalator introduced in [71] is directly related to the nearby canard orbits and so-called ribbons of the attracting slow manifolds. Such ribbons are separate surfaces, bounded by pairs of canard orbits that exhibit the same number of SAOs. We use the term twin canard orbits to describe a pair of canard orbits that bound a ribbon. We find a very similar structure of ribbons and associated twin canard orbits in the four-dimensional Hodgkin-Huxley model (32).
Similar to the approach taken in [38], we extend S a ε by the flow up to the threedimensional cross section Σ 1 , which was chosen as part of the selected boundary conditions for computing S s ε and associated canard orbits; see equations (45), (46), and (48). Figure 18 illustrates a selection of ribbons of the attracting slow manifold S a ε for system (32). Panel (a) shows five ribbons R 4 -R 8 of the extended attracting slow manifold S a ε in five different colors. The inset (a2) is an enlargement near the line L a 2 where the ribbons are easily distinguished. Each ribbon R i is computed individually as a family of orbit segments that start from L a 2 , make i SAOs, and end at the cross section Σ 1 . Figure 18(b) shows ribbon R 6 (red surface) and the bounding canard orbits ξ 6 and ξ 6 (thick red curves). The inset (b2) is an enlargement near L a 2 . This ribbon R 6 is representative of the other ribbons R i , each of which is bounded by a pair of canard orbits ξ i and ξ i that exhibit i SAOs. We again refer to the canard orbits ξ i and ξ i as twin canard orbits. Moreover, each ribbon R i rotates as a surface around the fold curve F 1 and makes i SAOs before reaching Σ 1 . Notice also from Fig. 18 that each ribbon R i folds over sharply in the process.
As mentioned before in Sect. 5.1, we find more than 26 periodic orbits for I = 9.74 with different MMO signatures. We find that ribbons of the extended attracting slow manifold and associated twin canard orbits give great insight into the mechanism of these different MMO signatures. Figure 19 shows three different periodic MMOs with different signatures superimposed on the associated ribbons of the extended attracting slow manifold for (32) with I = 9.74. Panels (a), (b) and (c) show the stable periodic MMO with signature 1 6 (black orbit), and unstable periodic MMOs 1 7 and 1 8 (gray orbits), together with ribbons R 6 , R 7 and R 8 (colored surfaces), respectively. The insets show enlargements near the line L a 2 where the twin canard orbits ξ i and  Panel (b) shows the ribbon R 6 with the bounding twin canard orbits ξ 6 and ξ 6 (thick curves); the inset is again an enlargement ξ i (thick curves) are labeled. Starting from the lower part of a given periodic MMO 1 i in all panels, the MMO stays very close to the ribbon R i in between the bounding twin canard orbits ξ i and ξ i and makes i SAOs before making a large excursion back to its starting point. This mechanism of MMOs was also found for the autocatalator system in [38].

Continuation of Canard Orbits
We now investigate the evolution of twin canard orbits as ε is varied. Once a canard orbit ξ i has been found as a solution of the well-posed 2PBVP with boundary conditions (48), it can be continued in system parameters. Figure 20 shows the continuation of canard orbits ξ 0 -ξ 15 in ε. We start from the canard orbits that we find for ε = 0.0083 and continue them in the direction of decreasing ε. Panel (a) shows the L 2 -norm of the canard orbits versus ε. In slight abuse of notation, we use ξ i to denote the ε-dependent branches that correspond to canard orbits with i SAOs. We find that the branches of canard orbits ξ 0 -ξ 9 can all be continued to the singular limit, and the respective canard orbits all converge to the strong singular canard ξ s in the limit Each ribbon R i is bounded by its respective pair of twin canard orbits ξ i and ξ i as ε → 0, as predicted by the theory [34]. Because of this accumulation onto ξ s , the continuation of canard orbits with larger numbers of SAOs becomes a very delicate task numerically as ε decreases, which is why ξ 10 -ξ 15 do not reach the singular limit ε = 0. Panels (b)-(e) of Fig. 20 demonstrate the evolution of the canard orbit ξ 6 in projection onto the (h, v)-plane for ε = 0.0083, ε = 0.003, ε = 3.0 × 10 −5 and ε = 1.0 × 10 −8 , respectively. As the singular limit ε = 0 is approached, the mechanism for SAOs of the canard orbits is clearly directly related to the folded node p, as is expected from the theory [33,34]. More precisely, the SAOs occur near p and away from q. Also, the amplitudes of the SAOs become infinitely small as ε → 0. This is another indication that our computational technique is reliable for computing canard orbits and saddle slow manifolds. The transition of canard orbit ξ 6 in panels (b)-(e) is representative of the other canard orbits ξ 1 -ξ 15 . Figure 21(a) shows the continuation of canard orbits ξ 1 -ξ 16 ; the branches of these canard orbits all go toward ε = 0, but the orbit segments do not converge to the strong singular canard ξ s . This is illustrated in Fig. 21(b)-(e), where we demonstrate the evolution of the canard orbit ξ 6 as ε decreases to 0. Here, ξ 6 is shown in projection onto the (h, v)-plane for ε = 0.0083, ε = 0.003, ε = 3.0 × 10 −5 and ε = 3.52578 × 10 −8 , respectively. Panel (c) shows that the last SAO of the canard orbit ξ 6 has already become much larger than the other SAOs. In panel (d), the last oscillation of ξ 6 becomes even larger and involves concatenations of slow and fast segments. Note that the canard orbit ξ 6 now looks like a canard orbit ξ 5 followed by a fast return segment to a canard orbit ξ 0 . The amplitudes of the first five SAOs become infinitely small in the limit as ε → 0. As a result, in the limit ε = 0, the canard orbit ξ 6 is a concatenation of slow and fast segments of the reduced system (34) and the fast subsystem (33), respectively. This difference in evolution, as compared to ξ 6 , is representative for all twin canard orbits ξ 1 -ξ 16 , and was also found in the reduced three-dimensional Hodgkin-Huxley model; see [36]. Figure 22 shows the continuation of canard orbits ξ 16 -ξ 21 , ξ 36 , ξ 75 and ξ 111 , as well as ξ 17 -ξ 21 , ξ 36 , ξ 75 and ξ 111 . Panels (a) and (b) shows the L 2 -norm of the canard orbits versus ε. We find that the shown branches do not converge to the singular limit, as expected from the theory [34]. Instead, the L 2 -norm of each branch increases asymptotically at a finite value of ε. This also confirms that canard orbits ξ 10 − ξ 15 and ξ 10 − ξ 16 are qualitatively different and should indeed go to the singular limit of ε. We remark that the limiting values of ε, for a given pair ξ i and ξ i in Fig. 22(a)-(b), are close but different.
To understand this transition, panels (c1) and (d1) of Fig. 22 show the canard orbits ξ 111 and ξ 111 , respectively, in projection onto the (h, v)-plane for ε = 0.0083. Notice that the mechanism for SAOs of ξ 111 and ξ 111 is directly related to the saddle focus q rather than the folded node p. More specifically, these canard orbits approach q along the weak direction of the two-dimensional stable manifold W s (q) and make 111 SAOs around q, due to the rotational nature of the two-dimensional unstable manifold W u (q). The difference is that ξ 111 has one final larger oscillation after the passage near q; compare (c1) and (d1). As ε decreases during the continuation, both canard orbits approach W s (q) much more closely; see Fig. 22(c2) and (d2). As a result, the respective starting points of the canard orbits lie on L a at h ≈ −21.897, that is, far to the left, outside the range of the panels. As ε is decreased further, the twin canard orbits ξ 111 and ξ 111 converge to W s (q) and connect to q along the weak stable eigendirection in the limit. Nevertheless, past q, the canard orbit ξ 111 still features the final larger oscillation.
When the twin canard orbits ξ i and ξ i are continued for increasing ε, we find that the respective branches have folds but do not merge with one another. This is in contrast to the autocatalator system [38], where twin canard orbits are created at fold bifurcations. Further investigation of the difference in the nature of twin canard orbits is beyond the scope of this paper and left for future research.

Conclusions and Directions for Future Work
We introduced a general method for computing two-dimensional saddle slow manifolds and associated canard orbits in systems with two fast and two slow variables. Our approach is based on the idea of choosing boundary conditions that define submanifolds transverse to stable and unstable normal directions of the saddle slow manifold. Indeed, its novelty is to define and solve two-point boundary-value problems (2PBVPs) with boundary conditions corresponding to codimension-i and codimension-j hypersurfaces for any i and j , such that i + j = 4. We employed and tested the resulting computational methods for an extended normal form of a folded node. Taking into account the nature of the slow flow and its normal direction, we were able to compute different relevant parts of attracting and saddle slow manifolds in the neighborhood of a folded node. The flexibility of our method also enabled us to compute two-dimensional submanifolds of the three-dimensional stable and unstable manifolds of a saddle slow manifold.
Furthermore, we illustrated how the interaction between two-dimensional slow manifolds in R 4 is very similar to that found in R 3 ; however, their intersections are no longer structurally stable. This means that canard orbits can no longer be described as intersections of two-dimensional slow manifolds. Rather, canard orbits exist as trajectories that follow attracting and saddle slow manifolds for O(1) time on the slow time scale before exiting via the unstable manifolds of the saddle slow manifold. We introduced a computational approach for detecting canard orbits in R 4 by defining suitable 2PBVPs and employing a homotopy step to find a first solution. With this approach, we systematically detected canard orbits and continued them in the timescale parameter ε, as well as the eigenvalue ratio μ of the folded node. Our findings demonstrate that the canard orbits of the extended normal form in R 4 behave according to what the theory predicts for the three-dimensional normal form.
We demonstrated that our computational approach is also practical for studying complex slow-fast dynamics that inherently appear in applications and especially in neuroscience. Specifically, we investigated the dynamics of the full four-dimensional Hodgkin-Huxley model without applying any reduction. The bifurcation diagram of this system features various signatures of MMOs, which are organized by slow manifolds and associated canard orbits. With our approach, we were also able to compute relevant submanifolds of the attracting and saddle slow manifolds, as well as associated canard orbits in the four-dimensional phase space of the system. We showed that, for ε sufficiently large, the extended attracting slow manifold is subdivided into ribbons that are bounded by twin canard orbits in a similar structure to that found for the three-dimensional autocatalator model [38]. Continuation analysis of canard orbits in ε showed again that the relevant subset of these canard orbits converges to the strong singular canard as the singular limit is approached.
When continuing canard orbits in ε in the four-dimensional Hodgkin-Huxley model, we found that the twin canard orbit ξ i+1 , for 0 ≤ i < 16, has the same limit as the canard orbit ξ i , but with an additional large oscillation. When ξ i+1 can be continued all the way to the singular limit ε = 0, this large oscillation converges to the primary strong canard ξ s composed with a fast return segment. This convergence process has also been found for a selection of canard orbits in the three-dimensional Hodgkin-Huxley model [36]. On the other hand, canard orbits ξ i and ξ i+1 , for i ≥ 16, do not converge to the singular limit in the four-dimensional Hodgkin-Huxley model, but connect to the saddle focus of the system as ε is decreased.
Overall, the work presented in this paper shows that it is possible to provide comprehensive insight into the organization of MMOs by examining the geometry of two-dimensional slow manifolds, even when the number of fast variables exceeds one. The basis of this achievement lies in a versatile construction of suitable selections of 2PBVPs, adapted to the respective context. In this way, the organization of the phase space by slow manifolds of four-dimensional slow-fast systems can be studied efficiently and comprehensively Interactions between slow and invariant manifolds are known to play an important role in the organization of MMOs. In particular, intersections of repelling slow manifolds and unstable manifolds of saddle-focus equilibria may organize the large excursions of MMOs [72][73][74]. Such connecting orbits can be found with a homotopy technique similar to that used for detecting canard orbits in this paper. Understanding the nature of such intersections between different types of manifolds in the Hodgkin-Huxley model and other slow-fast systems in R 4 is an interesting direction for future work.
Many examples in applications involve more than two time scales [39,44,75]. Furthermore, many slow-fast systems have no explicit time-scale separation and may contain various regions with different splittings of time scales [4,5,44,76,77]. In such systems, it will be necessary to study different regions locally using a variety of 2PBVP setups. These regions then need to be connected globally in order to understand the recurrent dynamics. The computational methods presented in this paper may be of use also for the investigation of such more general systems.
for x ∈ {m, h, n}, where p m = 120, p h = 1, p n = 1. The other parameters are given in Table 1; see Sect. 5. Note that the value of p m is different from that given in [42] so that it agrees with the nondimensionalization performed in [1] and corresponds to the classical parameter values of the original Hodgkin-Huxley model [8].

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.