#### Abstract

In a previous paper we have proposed a new method for proving the existence of “canard solutions” for three- and four-dimensional singularly perturbed systems with only one* fast* variable which improves the methods used until now. The aim of this work is to extend this method to the case of four-dimensional singularly perturbed systems with two* slow* and two* fast* variables. This method enables stating a unique generic condition for the existence of “canard solutions” for such four-dimensional singularly perturbed systems which is based on the stability of* folded singularities* (*pseudo singular points* in this case) of the* normalized slow dynamics* deduced from a well-known property of linear algebra. This unique generic condition is identical to that provided in previous works. Application of this method to the famous coupled FitzHugh-Nagumo equations and to the Hodgkin-Huxley model enables showing the existence of “canard solutions” in such systems.

#### 1. Introduction

In the beginning of the eighties, Benoît and Lobry [1], Benoît [2], and then Benoît [3] in his PhD- thesis studied canard solutions in . In the article entitled “Systèmes Lents-Rapides dans et Leurs Canards,” Benoît [2, p. 170] proved the existence of canards solution for three-dimensional singularly perturbed systems with two* slow* variables and one* fast* variable while using “nonstandard analysis” according to a theorem which stated that canard solutions exist in such systems provided that the* pseudo singular point* (this concept has been originally introduced by Argémi [4]; see Section 2.8) of the* slow dynamics*, that is, of the* reduced vector field,* is of* saddle-*type. Nearly twenty years later, Szmolyan and Wechselberger [5] extended “Geometric Singular Perturbation Theory (see Fenichel [6, 7], O’Malley [8], Jones [9], and Kaper [10])” to canards problems in and provided a “standard version” of Benoît’s theorem [2]. Very recently, Wechselberger [11] generalized this theorem for -dimensional singularly perturbed systems with * slow* variables and * fast* (1). The methods used by Szmolyan and Wechselberger [5] and Wechselberger [11] require implementing a “desingularization procedure” which can be summarized as follows: first, they compute the* normal form* of such singularly perturbed systems which is expressed according to some coefficients ( and for dimension three and , , and for dimension four) depending on the functions defining the original vector field and their partial derivatives with respect to the variables. Secondly, they project the “desingularized vector field” (originally called “normalized slow dynamics” by Benoît [2, p. 166]) of such a* normal form* on the tangent bundle of the critical manifold. Finally, they evaluate the Jacobian of the projection of this “desingularized vector field” at the* folded singularity* (originally called* pseudo singular points* by Argémi [4, p. 336]). This leads Szmolyan and Wechselberger [5, p. 427] and Wechselberger [11, p. 3298] to a “classification of* folded singularities* (*pseudo singular points*).” Thus, they showed that for three-dimensional singularly perturbed systems such* folded singularity* is of* saddl-type* if the following condition is satisfied, while for four-dimensional singularly perturbed systems such* folded singularity* is of* saddle-type* if . Then, Szmolyan and Wechselberger [5, p. 439] and Wechselberger [11, p. 3304] established their Theorem . which states that “In the folded saddle and in the folded node case singular canards perturb to maximal canard for sufficiently small ” However, in their works neither Szmolyan and Wechselberger [5] nor Wechselberger [11] did not provide (to our knowledge) the expression of these constants ( and ) which are necessary to state the existence of canard solutions in such systems.

In a previous paper entitled “Canards Existence in Memristor’s Circuits” (see Ginoux and Llibre [12]) we first provided the expression of these constants and then showed that they can be directly determined starting from the* normalized slow dynamics* and not from the projection of the “desingularized vector field” of the* normal form*. This method enabled stating a unique “generic” condition for the existence of “canard solutions” for such three- and four-dimensional singularly perturbed systems which is based on the stability of* folded singularities* of the* normalized slow dynamics* deduced from a well-known property of linear algebra. This unique condition which is completely identical to that provided by Benoît [2] and then by Szmolyan and Wechselberger [5] and finally by Wechselberger [11] is “generic” since it is exactly the same for singularly perturbed systems of dimensions three and four with only one* fast* variable.

The aim of this work is to extend this method to the case of four-dimensional singularly perturbed systems with * slow* and * fast* variables. Since the dimension of the system is , such problem is known as “canards existence in ” Moreover, in this particular case, where , the* folded singularities* of Wechselberger [11, p. 3298] are nothing else but the* pseudo singular points* of the late Argémi [4] as we will see below. Following the previous works, we show that for such four-dimensional singularly perturbed systems* pseudo singular points* are of* saddle-type* if . Then, according to Theorem . of Wechselberger [11, p. 3304] we provide the expression of this constant which is necessary to establish the existence of canard solutions in such systems. So, we can state that the condition for existence of canards in such is “generic” since it is exactly the same for singularly perturbed systems of dimensions three and four with only one* fast* variable.

The outline of this paper is as follows. In Section 1, definitions of singularly perturbed system, critical manifold, reduced system, “constrained system,” canard cycles, folded singularities, and pseudo singular points are recalled. The method proposed in this paper is presented in Section 2 for the case of four-dimensional singularly perturbed systems with two* fast* variables. In Section 3, applications of this method to the famous coupled FitzHugh-Nagumo equations and to the Hodgkin-Huxley model enable showing the existence of “canard solutions” in such systems.

#### 2. Definitions

##### 2.1. Singularly Perturbed Systems

According to Tikhonov [13], Jones [9], and Kaper [10]* singularly perturbed systems* are defined aswhere , , , and the prime denotes differentiation with respect to the independent variable . The functions and are assumed to be functions (in certain applications these functions will be supposed to be , ) of , , and in , where is an open subset of and is an open interval containing .

In the case when , that is, is a small positive number, the variable is called* slow* variable and is called* fast* variable. Using Landau’s notation, represents a function of and such that is bounded for positive going to zero, uniformly for in the given domain.

In general we consider that evolves at an rate while evolves at an * slow* rate. Reformulating system (1) in terms of the rescaled variable , we obtain

The dot represents the derivative with respect to the new independent variable .

The independent variables and are referred to the* fast* and* slow* times, respectively, and (1) and (2) are called the* fast* and* slow* systems, respectively. These systems are equivalent whenever , and they are labeled* singular perturbation problems* when . The label “singular” stems in part from the discontinuous limiting behavior in system (1) as .

##### 2.2. Reduced Slow System

In such case system (2) leads to a differential-algebraic system (DAE) called* reduced slow system* whose dimension decreases from to . Then, the* slow* variable partially evolves in the submanifold called the* critical manifold* (it represents the approximation of the slow invariant manifold, with an error of ). The* reduced slow system* is

##### 2.3. Slow Invariant Manifold

The* critical manifold* is defined by

Such a normally hyperbolic invariant manifold (4) of the* reduced slow system* (3) persists as a locally invariant* slow manifold* of the full problem (1) for sufficiently small. This locally* slow invariant manifold* is close to the* critical manifold*.

When is invertible, thanks to the Implicit Function Theorem, is given by the graph of a function for , where is a compact, simply connected domain and the boundary of is a -dimensional submanifold (the set is overflowing invariant with respect to (2) when ; see Kaper [10] and Jones [9]).

According to Fenichel [6, 7] theory if is sufficiently small, then there exists a function defined on such that the manifold,is locally invariant under the flow of system (1). Moreover, there exist perturbed local stable (or attracting) and unstable (or repelling) branches of the* slow invariant manifold *. Thus, normal hyperbolicity of is lost via a saddle-node bifurcation of the* reduced slow system* (3). Then, it gives rise to solutions of “canard” type.

##### 2.4. Canards: Singular Canards and Maximal Canards

A* canard* is a solution of singularly perturbed dynamical system (1) following the* attracting* branch of the* slow invariant manifold*, passing near a bifurcation point located on the fold of this* slow invariant manifold*, and then following the* repelling* branch of the* slow invariant manifold*.

A* singular canard* is a solution of* reduced slow system* (3) following the* attracting* branch of the* critical manifold*, passing near a bifurcation point located on the fold of this* critical manifold*, and then following the* repelling* branch of the* critical manifold*.

A* maximal canard* corresponds to the intersection of the attracting and repelling branches of the slow manifold in the vicinity of a nonhyperbolic point.

According to Wechselberger [11, p. 3302],

“Such a maximal canard defines a family of canards nearby which are exponentially close to the maximal canard, i.e., a family of solutions of (1) that follow an attracting branch of the slow manifold and then follow, rather surprisingly, a repelling/saddle branch of the slow manifold for a considerable amount of slow time. The existence of this family of canards is a consequence of the non-uniqueness of and . However, in the singular limit , such a family of canards is represented by a unique singular canard.”

Canards are special class of solutions of singularly perturbed dynamical systems for which normal hyperbolicity is lost. Canards in singularly perturbed systems with two or more slow variables , and one fast variable , are robust, since maximal canards generically persist under small parameter changes (see Benoît [2, 14], Szmolyan and Wechselberger [5], and Wechselberger [11, 15]).

##### 2.5. Constrained System

In order to characterize the “slow dynamics,” that is, the slow trajectory of the* reduced slow system* (3) (obtained by setting in (2)), Takens [16] introduced the “constrained system” defined as follows:

Since, according to Fenichel [6, 7], the* critical manifold * may be considered as locally invariant under the flow of system (1), we have

By replacing by leads to

This justifies the introduction of the* constrained system*.

Now, let denote the adjoint of the matrix which is the transpose of the cofactor matrix ; then while multiplying the left-hand side of (6) by the inverse matrix obtained by the adjoint method we have

##### 2.6. Normalized Slow Dynamics

Then, by rescaling the time by setting we obtain the following system which has been called by Benoît [2, p. 166] “normalized slow dynamics”:where the overdot now denotes the time derivation with respect to .

Let us notice that Argémi [4] proposed to rescale time by setting in order to keep the same flow direction in (10) as in (9).

##### 2.7. Desingularized Vector Field

By application of the Implicit Function Theorem, let suppose that we can explicitly express from (4), say, without loss of generality, as a function of the other variables. This implies that is locally the graph of a function over the base , where . Thus, we can span the “normalized slow dynamics” on the tangent bundle at the* critical manifold * at the* pseudo singular point*. This leads to the so-called* desingularized vector field*:

##### 2.8. Pseudo Singular Points and Folded Singularities

As recalled by Guckenheimer and Haiduc [17, p. 91],* pseudo singular points* have been introduced by the late Argémi [4] for low-dimensional singularly perturbed systems and are defined as singular points of the “normalized slow dynamics” (10). Twenty-three years later, Szmolyan and Wechselberger [5, p. 428] called such* pseudo singular points folded singularities*. In a recent publication entitled “A Propos de Canards” Wechselberger [11, p. 3295] proposed to define such singularities for -dimensional singularly perturbed systems with * slow* variables and * fast* as the solutions of the following system:

Thus, for dimensions higher than three, his concept encompasses that of Argémi. Moreover, Wechselberger [11, p. 3296] proved that* folded singularities* form a -dimensional manifold. Thus, for the* folded singularities* are nothing else than the* pseudo singular points* defined by Argémi [4], while for the* folded singularities* are no more points but a -dimensional manifold. Moreover, let us notice on one hand that the original system (1) includes variables and on the other hand that system (12) comprises equations. However, in the particular case , two equations of system (12) are linearly dependent. So, such system only comprises equations. So, all the variables (unknowns) of system (12) can be determined. The solutions of this system are called* pseudo singular points*. We will see in Section 3 that the stability analysis of these* pseudo singular points* will give rise to a condition for the existence of canard solutions in the original system (1).

#### 3. Four-Dimensional Singularly Perturbed Systems with Two Fast Variables

Four-dimensional* singularly perturbed dynamical system* (2) with * slow* variables and * fast* may be written aswhere , , , and the functions and are assumed to be functions of .

##### 3.1. Critical Manifold

The critical manifold equation of system (13a), (13b), (13c), and (13d) is defined by setting in ((13c) and (13d)). Thus, we obtain

By application of the Implicit Function Theorem, let us suppose that we can explicitly express from ((14a) and (14b)), say, without loss of generality, and as functions of the other variables:

##### 3.2. Constrained System

The* constrained system* is obtained by equating with zero the time derivative of :

Equations (16a) and (16b) may be written as

By solving the system of (17a) and (17b) with two unknowns we find

So, we have the following constrained system:

##### 3.3. Normalized Slow Dynamics

By rescaling the time by setting we obtain the “normalized slow dynamics”:where the overdot now denotes the time derivation with respect to .

##### 3.4. Desingularized System on the Critical Manifold

Then, since we have supposed that and may be explicitly expressed as functions of the other variables ((15a) and (15b)), they can be used to project the normalized slow dynamics (20) on the tangent bundle of the critical manifold. So, we have

##### 3.5. Pseudo Singular Points

*Pseudo singular points* are defined as singular points of the “normalized slow dynamics,” that is, as the set of points for which we have

*Remark 1. *
Let us notice on one hand that (22b) and (22c) are linearly dependent and on the other hand that contrary to previous works we do not use the “desingularized vector field” (21) but the “normalized slow dynamics” (20).

The Jacobian matrix of system (20) reads

##### 3.6. Extension of Benoît’s Generic Hypothesis

Without loss of generality, it seems reasonable to extend Benoît’s generic hypotheses introduced for the three-dimensional case to the four-dimensional case. So, first, let us suppose that by a “standard translation” the* pseudo singular point* can be shifted at the origin and that by a “standard rotation” of -axis that the* slow manifold* is tangent to ()-hyperplane, so we have

Then, let us make the following assumptions for the nondegeneracy of the* folded singularity*:

According to these generic hypotheses ((24)-(25)), the Jacobian matrix (23) readswhere

Thus, we have the following Cayley-Hamilton eigenpolynomial associated with such Jacobian matrix (26) evaluated at the* pseudo singular point*, that is, at the origin:where is the sum of all first-order diagonal minors of , that is, the trace of the Jacobian matrix , represents the sum of all second-order diagonal minors of , and represents the sum of all third-order diagonal minors of . It appears that since one row of the Jacobian matrix (26) is null. So, the eigenpolynomial reduces to

But, according to Wechselberger [11], vanishes at a* pseudo singular point* as it is easy to prove it. So, the eigenpolynomial (29) is reduced to

Let be the eigenvalues of eigenpolynomial (30) and let us denote by the obvious double root of this polynomial. We havewhere is is the sum of all first-order diagonal minors of , that is, the trace of the Jacobian matrix , and represents the sum of all second-order diagonal minors of . Thus, the* pseudo singular point* is of saddle-type* iff* the following conditions and are verified:

Condition is systematically satisfied provided that condition is verified. Thus, the* pseudo singular point* is of saddle-type* iff *.

##### 3.7. Canard Existence in

Following the works of Wechselberger [11] it can be stated, while using a standard polynomial change of variables, that any -dimensional singularly perturbed systems with * slow* variables () and * fast* () (1) can be transformed into the following “normal form”:

We establish in the Appendix for any four-dimensional singularly perturbed systems (13a), (13b), (13c), and (13d) with * slow* and * fast* variables that

Thus, in his paper Wechselberger [11, p. 3304] provided in the framework of “standard analysis” a generalization of Benoît’s theorem [2] for any -dimensional singularly perturbed systems with * slow* variables () and * fast* (). According to his theorem presented below as Theorem 2 he proved the existence of canard solutions for the original system (1).

Theorem 2.
*In the folded saddle case of system (33) singular canards perturb to maximal canards solutions for sufficiently small .*

*Proof. *
See Wechselberger [11].

Since our method does not use the “desingularized vector field” (21) but the “normalized slow dynamics” (20), we have the following proposition.

Proposition 3.
*
If the normalized slow dynamics (20) has a pseudo singular point of saddle-type, that is, if the sum of all second-order diagonal minors of the Jacobian matrix of the normalized slow dynamics (20) evaluated at the pseudo singular point is negative, that is, if then, according to Theorem 2, system (13a), (13b), (13c), and (13d) exhibits a canard solution which evolves from the attractive part of the slow manifold towards its repelling part.*

*Proof. *
By making some smooth changes of time and smooth changes of coordinates (see the Appendix) we brought system (13a), (13b), (13c), and (13d) to the following “normal form”:Then, we deduce that the condition for the* pseudo singular point* to be of saddle-type is . According to (32) it is easy to verify thatSo, the condition for which the* pseudo singular point* is of saddle-type, that is, , is identical to that proposed by Wechselberger [11, p. 3298] in his theorem; that is, .

So, Proposition 3 can be used to state the existence of canard solution for such systems. Application of Proposition 3 to the coupled FitzHugh-Nagumo equations, presented in Section 4, which is a four-dimensional singularly perturbed system with two* slow* and two* fast* variables will enable proving, as many previous works such as those of Tchizawa and Campbell [18] and Tchizawa [19–24], the existence of “canard solutions” in such system. According to Tchizawa [25], it is very important to notice on one hand that the fast equation has 2 dimensions in the system and on the other hand that the fast system can give attractive, repulsive, or attractive-repulsive part at each* pseudo singular point*. Then, Tchizawa [25] has established that the jumping direction can be shown using the eigenvectors. In the same way we will find again the results of Rubin and Wechselberger [26] concerning the existence of “canard solutions” in the Hodgkin-Huxley model but with a set of more realistic parameters used in Chua et al. [27, 28].

#### 4. Coupled FitzHugh-Nagumo Equations

The FitzHugh-Nagumo model [29, 30] is a simplified version of the Hodgkin-Huxley model [31] which models in a detailed manner activation and deactivation dynamics of a spiking neuron. By coupling two FitzHugh-Nagumo models Tchizawa and Campbell [18] and Tchizawa [19, 24] obtained the following four-dimensional singularly perturbed system with two* slow* and two* fast* variables:where and is the “canard parameter” or “duck parameter” while is a scale factor.

##### 4.1. Slow Manifold and Constrained System

The slow manifold equation of system (37a), (37b), (37c), and (37d) is defined by setting in (37c) and (37d). Thus, we obtain

##### 4.2. Normalized Slow Dynamics

Then, by rescaling the time by setting we obtain the “normalized slow dynamics”:

##### 4.3. Pseudo Singular Points

From (22a), (22b), (22c), (22d), and (22e), the* pseudo singular points* of system (37a), (37b), (37c), and (37d) are defined by

According to Tchizawa and Campbell [18] and Tchizawa [19, 20], there are six* pseudo singular points* with the last four depending on the parameter :

##### 4.4. Canard Existence in Coupled FitzHugh-Nagumo Equations

The Jacobian matrix of system (39) evaluated at the* pseudo singular points* (41a) reads

*Remark 4. *
Although the* pseudo singular points* have not been shifted at the origin extension of Benoît’s generic hypotheses (24)-(25) are satisfied. In other words, we have .

According to (31) we find that

Thus, according to Proposition 3, the* pseudo singular points* are of saddle-type if and only if:

So, we have conditions and :

Let us choose arbitrarily as the “canard parameter” or “duck parameter.” Obviously, it appears that the condition is still satisfied. Finally, the* pseudo singular points* are of saddle-type if and only if we have

*Remark 5. *
Let us notice that the* pseudo singular points* are of node-type if as stated by Tchizawa and Campbell [18] and Tchizawa [19, 20].

The Jacobian matrix of system (39) evaluated at the* pseudo singular points* (41b) reads

*Remark 6. *
Although the* pseudo singular points* have not been shifted at the origin extension of Benoît’s generic hypotheses (24)-(25) are satisfied. In other words, we have .

According to (31) we find that

Thus, according to Proposition 3, the* pseudo singular points* are of saddle-type if and only if

So, we have conditions and :

Let us choose arbitrarily as the “canard parameter” or “duck parameter.” Obviously, it appears that if the condition is verified, then the condition is* de facto* satisfied. Finally, the* pseudo singular points* are of saddle-type if and only if we have

*Remark 7. *
Because of the symmetry of these coupled FitzHugh-Nagumo equations, the Jacobian matrix of system (39) evaluated at the* pseudo singular points* (41c) provides the same result as just above.

#### 5. Hodgkin-Huxley Model

The original Hodgkin-Huxley model [31] is described by the following system of four nonlinear ordinary differential equations:

where:

The first equation (52a) results from the application of Kirchhoff’s law to the space clamped squid giant axon. Thus, the total membrane current for which represents the specific membrane capacity and the displacement of the membrane potential from its resting value is equal to the sum of the following intrinsic currents:where is a delayed rectifier potassium current, is fast sodium current, and is the “leakage current.” The parameter is the total membrane current density, inward positive, that is, the total current injected into the space clamped squid giant axon, and , , and are the equilibrium potentials of potassium, sodium, and “leakage current,” respectively. The maximal specific conductance degrees of the ionic currents are denoted by , , and , respectively. Functions and are gates’ opening and closing rates depending on . Variable denotes the activation of the sodium current, variable the inactivation of the sodium current, and variable the activation of the potassium current. These dimensionless gating variables vary in the range .

Let us notice that the variables and symbols in ((52a), (52b), (52c), and (52d) and (53a), (53b), (53c), (53d), (53e), and (53f)) originally chosen by Hodgkin and Huxley and are different from those found in recent literatures, where the reference polarity of the voltage and the reference direction of the current are defined as the negative of the voltages and currents. We have opted to adopt the reference assumption in Hodgkin and Huxley [31] for ease in comparison of our results with those from Hodgkin and Huxley (for more details see Chua et al. [27, 28]). The parameter values are exactly those chosen in the original Hodgkin-Huxley [31] works:

According to Suckley and Biktashev [32] and Suckley [33], dimensionless functions , , and called gates’ instant equilibrium values, that is, steady-state relation for gating variables , , and , respectively, as well as , , and called gates dynamics timescales in , that is, time constant for gating variables , , and , respectively, may be defined as follows:

By using (56a), (56b), (56c), (56d), (56e), and (56f), the original Hodgkin-Huxley model [31] reads

Now, in order to apply the* singular perturbation method* to the Hodgkin-Huxley model, two small multiplicative parameters are introduced. According to Suckley and Biktashev [32], Suckley [33], and Rubin and Wechselberger [26], the existence of two different timescales of evolution for couples of dynamic variables () and () enables justifying such an introduction. So, in order to differentiate* slow* variables from* fast* variables, Suckley and Biktashev [32], Suckley [33], and Rubin and Wechselberger [26] have plotted the inverse of “time constant for gating variable ,” that is, according to with . In Figure 1, they have been plotted for the original functions and (53a). However, let us notice that this plot is exactly the same as those presented by Rubin and Wechselberger [26] (Figure ) for a nondimensionalized three-dimensional Hodgkin-Huxley singularly perturbed system obtained after the following variable changes: and ; then and finally .

Figure 1 shows a plot of the functions according to with over the physiological range. We observe that is of an order of magnitude bigger than and , which are of comparable size. Indeed, we can deduce that the values of times scales are approximately while . Then, it appears that corresponds to the* fast* variable while and correspond to* slow* variables. Moreover, since the activation of the sodium channel is directly related to the dynamics of the membrane (action) potential , Rubin and Wechselberger [26] consider that and evolve on the same* fast* timescale. So, the Hodgkin-Huxley model may be transformed into a* singularly perturbed system* with two timescales in which the* slow* variables are and the* fast* variables are .

So, according to Awiszus et al. [34], Suckley and Biktashev [32], Suckley [33], and Rubin and Wechselberger [26] small multiplicative parameters in the original vector field of the Hodgkin-Huxley equations (57a), (57b), (57c), and (57d) may be identified while factorizing the right-hand side of (57a) by and set:Other parameters are kept as for the original Hodgkin-Huxley model [31]:

Then, by posing , , and () () to be consistent with the notations of Section 3, we obtainwhere () () and .

Let us notice that the multiplicative parameter has been introduced artificially in (60c). This is due to the fact that it has been stated above that the timescale of variable , that is, , is tenth times greater than the timescale of variables and , that is, of variables and . Moreover, this parameter is identical to those used in (60d) since it has been also considered that and , that is, and , evolve on the same* fast* timescale.

According to the* Geometric Singular Perturbation Theory*, the zero-order approximation in of the* slow manifold* associated with the Hodgkin-Huxley model (60a), (60b), (60c), and (60d) is obtained by posing in (60c) and (60d). So, the* slow manifold* is given by

Then, the* fast foliation* is within the planes and .

The* fold curve* is defined as the location of the points, where , , and . For the Hodgkin-Huxley model (60a), (60b), (60c), and (60d), the* fold curve* is thus given by (61a) and (61b) and by the determinant of the Jacobian matrix of the following * fast foliation*:

The Jacobian matrix of the* fast foliation* (62a) and (62b) readswhere denotes the derivative with respect to . Then, taking into account (61b), that is, , we have

So, the determinant of the Jacobian matrix of the* fast foliation* (62a) and (62b) is

Thus, the condition for the* fold curve* is , which gives

Therefore