Canards Existence in FitzHugh-Nagumo and Hodgkin-Huxley Neuronal Models

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 to state 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. Applications of this method to the famous coupled FitzHugh-Nagumo equations and to the Hodgkin-Huxley model enables to show the existence of"canard solutions"in such systems.


INTRODUCTION
In the beginning of the eighties, Benoît and Lobry [5], Benoît [6] and then Benoît [7] in his PhD-thesis studied canard solutions in R 3 . In the article entitled "Systèmes lents-rapides dans R 3 et leurs canards," Benoît [6, p. 170] proved the existence of canards solution for three-dimensional singularly perturbed systems with two slow variables and one fast variable while using "Non-Standard Analysis"according to a theorem which stated that canard solutions exist in such systems provided that the pseudo singular point 1 of the slow dynamics, i.e., of the reduced vector field is of saddle type. Nearly twenty years later, Szmolyan and Wechselberger [25] extended "Geometric Singular Perturbation Theory 2 " to canards problems in R 3 and provided a "standard version" of Benoît's theorem [6]. Very recently, Wechselberger [39] generalized this theorem for n-dimensional singularly perturbed systems with k slow variables and m fast (Eq. (1)). The method used by Szmolyan and Wechselberger [25] and Wechselberger [39] require to implement 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 (a and b for dimension three andã,b andc j 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 Eric Benoît [6, 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 José Argémi [1, p. 336]). This lead Szmolyan and Wechselberger [25, p. 427] and Wechselberger [39, 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 saddle type if the following condition is satisfied: a < 0 while for four-dimensional singularly perturbed systems such folded singularity is of saddle type ifã < 0. Then, Szmolyan and Wechselberger [25, p. 439] and Wechselberger [39, p. 3304] established their Theorem 4.1. which state 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 [25] nor Wechselberger [39] did not provide (to our knowledge) the expression of these constants (a 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 & Llibre [17]) 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 to state 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 [6] and then by Szmolyan and Wechselberger [25] and finally by Wechselberger [39] is "generic" since it is exactly the same for singularly perturbed systems of dimension three and four with only one fast variable.
The aim of this work is to extend this method to the case of fourdimensional singularly perturbed systems with k = 2 slow and m = 2 fast variables. Since the dimension of the system is m = k + m, such problem is known as "canards existence in R 2+2 ". Moreover, in this particular case where k = m = 2, the folded singularities of Wechselberger [39, p. 3298] are nothing else but the pseudo singular points of the late José Argémi [1] 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ã < 0. Then, according Theorem 4.1. of Wechselberger [39, 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ã < 0 for existence of canards in such R 2+2 is "generic" since it is exactly the same for singularly perturbed systems of dimension three and four with only one fast variable.
The outline of this paper is as follows. In Sec. 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 article is presented in Sec. 2 for the case of four-dimensional singularly perturbed systems with two fast variables. In Sec. 3, applications of this method to the famous coupled FitzHugh-Nagumo equations and to the Hodgkin-Huxley model enables to show the existence of "canard solutions" in such systems.
2. DEFINITIONS 2.1. Singularly perturbed systems According to Tikhonov [37], Jones [20] and Kaper [21] singularly perturbed systems are defined as: where x ∈ R k , y ∈ R m , ε ∈ R + , and the prime denotes differentiation with respect to the independent variable t ′ . The functions f and g are assumed to be C ∞ functions 3 of x, y and ε in U × I, where U is an open subset of R k × R m and I is an open interval containing ε = 0.
In the case when 0 < ε ≪ 1, i.e. ε is a small positive number, the variable x is called slow variable, and y is called fast variable. Using Landau's notation: O (ε p ) represents a function f of u and ε such that f (u, ε)/ε p is bounded for positive ε going to zero, uniformly for u in the given domain.
In general we consider that x evolves at an O (ε) rate; while y evolves at an O (1) slow rate. Reformulating system (1) in terms of the rescaled variable t = εt ′ , we obtain˙ The dot represents the derivative with respect to the new independent variable t.
The independent variables t ′ and t 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 ε = 0, and they are labeled singular perturbation problems when 0 < ε ≪ 1. The label "singular" stems in part from the discontinuous limiting behavior in system (1) as ε → 0.

Reduced slow system
In such case system (2) leads to a differential-algebraic system (D.A.E.) called reduced slow system whose dimension decreases from k + m = n to m. Then, the slow variable x ∈ R k partially evolves in the submanifold M 0 called the critical manifold 4 . The reduced slow system iṡ

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 O(ε) close to the critical manifold.
When D x f is invertible, thanks to the Implicit Function Theorem, M 0 is given by the graph of a C ∞ function x = G 0 ( y) for y ∈ D, where D ⊆ R k is a compact, simply connected domain and the boundary of D is According to Fenichel [12,15] theory if 0 < ε ≪ 1 is sufficiently small, then there exists a function G ( y, ε) defined on D such that the manifold is locally invariant under the flow of system (1). Moreover, there exist perturbed local stable (or attracting) M a and unstable (or repelling) M r branches of the slow invariant manifold M ε . Thus, normal hyperbolicity of M ε is lost via a saddle-node bifurcation of the reduced slow system (3). Then, it gives rise to solutions of "canard" type.

Canards, singular canards and maximal canards
A canard is a solution of a singularly perturbed dynamical system (1) following the attracting branch M a 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 M r of the slow invariant manifold.
A singular canard is a solution of a reduced slow system (3) following the attracting branch M a,0 of the critical manifold, passing near a bifurcation point located on the fold of this critical manifold, and then following the repelling branch M r,0 of the critical manifold.
A maximal canard corresponds to the intersection of the attracting and repelling branches M a,ε ∩ M r,ε of the slow manifold in the vicinity of a non-hyperbolic point.
According to Wechselberger [39, 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 M a,ε of the slow manifold and then follow, rather surprisingly, a repelling/saddle branch M r,ε 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 M a,ε and M r,ε . However, in the singular limit ε → 0, such a family of canards is represented by a unique singular canard." Canards are a special class of solution of singularly perturbed dynamical systems for which normal hyperbolicity is lost. Canards in singularly perturbed systems with two or more slow variables ( x ∈ R k , k 2) and one fast variable ( y ∈ R m , m = 1) are robust, since maximal canards generically persist under small parameter changes 6 .

Constrained system
In order to characterize the "slow dynamics", i.e. the slow trajectory of the reduced slow system (3) (obtained by setting ε = 0 in (2)), Floris Takens [28] introduced the "constrained system" defined as follows: Since, according to Fenichel [12,15], the critical manifold g ( x, y, 0) may be considered as locally invariant under the flow of system (1), we have: By replacing˙ x by f ( x, y, 0) leads to: This justifies the introduction of the constrained system. Now, let adj(D y g) denote the adjoint of the matrix D y g which is the transpose of the co-factor matrix D y g, then while multiplying the left hand side of (6) by the inverse matrix (D y g) −1 obtained by the adjoint method we have:˙ 2.6. Normalized slow dynamics Then, by rescaling the time by setting t = −det(D y g)τ we obtain the following system which has been called by Eric Benoît [6, p. 166] "normalized slow dynamics": where the overdot now denotes the time derivation with respect to τ . Let's notice that José Argémi [1] proposed to rescale time by setting t = −det(D y g)sgn(det(D y g))τ in order to keep the same flow direction in (8) as in (7).

Desingularized vector field
By application of the Implicit Function Theorem, let suppose that we can explicitly express from Eq. (4), say without loss of generality, x 1 as a function φ 1 of the other variables. This implies that M 0 is locally the graph of a function φ 1 : R k → R m over the base U = ( χ, y) where χ = (x 2 , x 3 , ..., x k ). Thus, we can span the "normalized slow dynamics" on the tangent bundle at the critical manifold M 0 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 [18, p. 91], pseudo-singular points have been introduced by the late José Argémi [1] for low-dimensional singularly perturbed systems and are defined as singular points of the "normalized slow dynamics" (8). Twenty-three years later, Szmolyan and Wechselberger [25, p. 428] called such pseudo singular points, folded singularities. In a recent publication entitled "A propos de canards" Wechselberger [39, p. 3295] proposed to define such singularities for n-dimensional singularly perturbed systems with k slow variables and m fast as the solutions of the following system: Thus, for dimensions higher than three, his concept encompasses that of Argémi. Moreover, Wechselberger [39, p. 3296] proved that folded singularities form a (k − 2)-dimensional manifold. Thus, for k = 2 the folded singularities are nothing else than the pseudo singular points defined by Argémi [1]. While for k 3 the folded singularities are no more points but a (k − 2)-dimensional manifold. Moreover, let's notice on the one hand that the original system (1) includes n = k + m variables and on the other hand, that the system (10) comprises p = 2m + 1 equations. However, in the particular case k = m = 2, two equations of the system (10) are linearly dependent. So, such system only comprisesp = 2m = 2k equations. So, all the variables (unknowns) of system (10) can be determined. The solutions of this system are called pseudo singular points. We will see in the next Sec. 2 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).

FOUR-DIMENSIONAL SINGULARLY PERTURBED SYSTEMS WITH TWO FAST VARIABLES
A four-dimensional singularly perturbed dynamical system (2) with k = 2 slow variables and m = 2 fast may be written as: and the functions f i and g i are assumed to be C 2 functions of (x 1 , x 2 , y 1 , y 2 ).

Critical Manifold
The critical manifold equation of system (11) is defined by setting ε = 0 in Eqs. (11c & 11d). Thus, we obtain: By application of the Implicit Function Theorem, let suppose that we can explicitly express from Eqs. (12a & 12b), say without loss of generality, x 1 and y 1 as functions of the others variables:

Constrained system
The constrained system is obtained by equating to zero the time derivative of g 1,2 (x 1 , x 2 , y 1 , y 2 ): Eqs. (14a & 14b) may be written as: By solving the system of two equations (15a & 15b) with two unknowns (ẏ 1 ,ẏ 2 ) we find: So, we have the following constrained system: 3.3. Normalized slow dynamics By rescaling the time by setting t = −det J (y1,y2) τ 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 x 1 and y 1 may be explicitly expressed as functions of the others variables (13a & 13b), they can be used to project the normalized slow dynamics (18) 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", i.e. as the set of points for which we have: Remark 1. Let's notice on the one hand that Eqs. (20b) & (20c) are linearly dependent and on the other hand that contrary to previous works we don't use the "desingularized vector field" (19) but the "normalized slow dynamics" (18).
The Jacobian matrix of system (18) reads:

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's suppose that by a "standard translation" the pseudo singular point can be shifted at the origin O(0, 0, 0, 0) and that by a "standard rotation" of y 1 -axis that the slow manifold is tangent to (x 2 , x 3 , y 1 )hyperplane, so we have Then, let's make the following assumptions for the non-degeneracy of the folded singularity: where Thus, we have the following Cayley-Hamilton eigenpolynomial associated with such a Jacobian matrix (24) evaluated at the pseudo singular point, i.e., at the origin: where σ 1 = T r(J) is the sum of all first-order diagonal minors of J, i.e., the the trace of the Jacobian matrix J, σ 2 represents the sum of all secondorder diagonal minors of J and σ 3 represents the sum of all third-order diagonal minors of J. It appears that σ 4 = |J| = 0 since one row of the Jacobian matrix (24) is null. So, the eigenpolynomial reduces to: But, according to Wechselberger [39], σ 3 vanishes at a pseudo singular point as it's easy to prove it. So, the eigenpolynomial (26) is reduced to Let λ i be the eigenvalues of the eigenpolynomial (27) and let's denote by λ 3,4 = 0 the obvious double root of this polynomial. We have: where σ 1 = T r(J (F1,F2,G1,G2) ) = p is is the sum of all first-order diagonal minors of J (F1,F2,G1,G2) , i.e. the trace of the Jacobian matrix J (F1,F2,G1,G2) and σ 2 = 3 i=1 J ii (F1,F2,G1,G2) = q represents the sum of all second-order diagonal minors of J (F1,F2,G1,G2) . Thus, the pseudo singular point is of saddle-type iff the following conditions C 1 and C 2 are verified: Condition C 1 is systematically satisfied provided that condition C 2 is verified. Thus, the pseudo singular point is of saddle-type iff q < 0.

Canard existence in R 2+2
Following the works of Wechselberger [39] it can be stated, while using a standard polynomial change of variables, that any n-dimensional singularly perturbed systems with k slow variables (k 2) and m fast (m 1) (1) can be transformed into the following "normal form": We establish in Appendix A for any four-dimensional singularly perturbed systems (11) with k = 2 slow and m = 2 fast variables that Thus, in his paper Wechselberger [39, p. 3304] provided in the framework of "standard analysis" a generalization of Benoît's theorem [6] for any ndimensional singularly perturbed systems with k slow variables (k 2) and m fast (m 1). According to his Theorem 4.1 presented below he proved the existence of canard solutions for the original system (1).

Theorem 2.
In the folded saddle case of system (30) singular canards perturb to maximal canards solutions for sufficiently small ε ≪ 1.
Since our method doesn't use the "desingularized vector field" (19) but the "normalized slow dynamics" (18), we have the following proposition: If the normalized slow dynamics (18) has a pseudo singular point of saddle type, i.e. if the sum σ 2 of all second-order diagonal minors of the Jacobian matrix of the normalized slow dynamics (18) evaluated at the pseudo singular point is negative, i.e. if σ 2 < 0 then, according to Theorem 2, system (11) 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 Appendix A) we brought the system (11) to the following "normal form": x 1 x 2 y 2 , Then, we deduce that the condition for the pseudo singular point to be of saddle type isã < 0. According to Eqs. (29) it is easy to verify that σ 1 = T r(J (F1,F2,G1,G2) ) = λ 1 + λ 2 = −bc, So, the condition for which the pseudo singular point is of saddle type, i.e. σ 2 < 0 is identical to that proposed by Wechselberger [39, p. 3298] in his theorem, i.e.ã < 0. So, Prop. 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 Sec. 4, which is a four-dimensional singularly perturbed system with two slow and two fast variables will enable to prove, as many previous works such as those of Tchizawa & Campbell [30] and Tchizawa [30,31,32,33,34,35], the existence of "canard solutions" in such system. According to Tchizawa [36], it is very important to notice, on the one hand that the fast equation has 2-dimensional in the system R 2+2 and, on the other hand that the fast system can give attractive, repulsive or attractive-repulsive at each pseudo singular point. Then, Tchizawa [36] has established that the jumping direction can be shown using the eigenvectors. In the same way we will find again the results of Rubin et al. [24] concerning the existence of "canard solutions" in the Hodgkin-Huxley model but with a set of more realistic parameters used in Chua et al. [10,11].

COUPLED FITZHUGH-NAGUMO EQUATIONS
The FitzHugh-Nagumo model [16,22] is a simplified version of the Hodgkin-Huxley model [19] which models in a detailed manner activation and deactivation dynamics of a spiking neuron. By coupling two FitzHugh-Nagumo models Tchizawa & Campbell [29] and Tchizawa [30,35] obtained the following four-dimensional singularly perturbed system with two slow and two fast variables: where 0 < ε ≪ 1 and b is the "canard parameter" or "duck parameter" while c is a scale factor.

Slow manifold and contrained system
The slow manifold equation of system (31) is defined by setting ε = 0 in Eqs. (31c & 31d). Thus, we obtain:

Canard existence in coupled FitzHugh-Nagumo equations
The Jacobian matrix of system (33) evaluated at the pseudo singular points (35a) reads: Remark 4. Although the pseudo singular points have not been shifted at the origin extension of Benoît's generic hypotheses (22)(23) are satisfied. In other words, we have σ 4 = σ 3 = 0.
According to Eqs. (28) we find that: Thus, according to Prop. 3, the pseudo singular points are of saddle-type if and only if: So, we have the following conditions C 1 and C 2 : Let's choose arbitrarily b as the "canard parameter" or "duck parameter". Obviously, it appears that the condition C 1 is still satisfied. Finally, the pseudo singular points are of saddle-type if and only if we have: Remark 5. Let's notice that the pseudo singular points are of nodetype if − 3 4 < b < 0 as stated by Tchizawa & Campbell [29] and Tchizawa [30,31].
According to Eqs. (28) we find that: Thus, according to Prop. 3, the pseudo singular points are of saddle-type if and only if: So, we have the following conditions C 1 and C 2 : Let's choose arbitrarily b as the "canard parameter" or "duck parameter". Obviously, it appears that if the condition C 2 is verified then the condition C 1 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 this coupled FitzHugh-Nagumo equations, the Jacobian matrix of system (33) evaluated at the pseudo singular points (35c) provides the same result as just above.

HODGKIN-HUXLEY MODEL
The original Hodgkin-Huxley model [19] is described by the following system of four nonlinear ordinary differential equations: where: The first equation (44a) results from the application of Kirchhoff's law to the space clamped squid giant axon. Thus, the total membrane current C M dV /dt for which C M represents the specific membrane capacity and V the displacement of the membrane potential from its resting value, is equal to the sum of the following intrinsic currents: where I K is a delayed rectifier potassium current, I N a is fast sodium current and I L is the "leakage current". The parameter I is the total membrane current density, inward positive, i.e. the total current injected into the space clamped squid giant axon and V K , V N a and V L are the equilibrium potentials of potassium, sodium and "leakage current" respectively. The maximal specific conductances of the ionic currents are denotedḡ K , g N a andḡ L respectively. Functions α n,m,h and β n,m,h are gates' opening and closing rates depending on V . Variable m denotes the activation of the sodium current, variable h the inactivation of the sodium current and variable n the activation of the potassium current. These dimensionless gating variables vary between [0, 1].
Let's notice that the variables and symbols in Eqs. According to Suckley and Biktashev [26] and Suckley [27], dimensionless functionsn,h andm called gates' instant equilibrium values, i.e., steadystate relation for gating variable n, h and m respectively as well as τ n , τ h and τ m called gates dynamics time scales in ms, i.e., time constant for gating variable n, h and m respectively may be defined as follows: By using Eqs. 46, the original Hodgkin-Huxley model [19] reads: Now, in order to apply the singular perturbation method to the Hodgkin-Huxley model, two small multiplicative parameters ε ≪ 1 are introduced. According to Suckley and Biktashev [26], Suckley [27] and Rubin and Wechselberger [24], the existence of two different time scales of evolution for the couples of dynamic variables (n, h) and (m, V ) enables to justify such an introduction. So, in order to differentiate slow variables from fast variables, Suckley and Biktashev [26], Suckley [27] and Rubin and Wechselberger [24] have plotted the inverse of "time constant for gating variable i", i.e., τ i −1 according to V with i = n, h, m. In Fig. 1, they have been plotted for the original functions α i and β i (Eqs. 45a). However, let's notice that this plot is exactly the same as those presented by Rubin and Wechselberger [24] (Fig. 1) for a nondimensionalized three-dimensional Hodgkin-Huxley singularly perturbed system obtained after the following variable changes: V → −V andĪ → −Ī, then V → V + 65 and finally V → V /100.   1 shows a plot of the functions τ i −1 according to V with i = n, h, m over the physiological range. We observe that τ m −1 is of an order of magnitude bigger than τ h −1 and τ n −1 , which are of comparable size. Indeed, we can deduce that the values of times scales are approximately Then, it appears that m corresponds to the fast variable while n and h correspond to slow variables. Moreover, since the activation of the sodium channel m is directly related to the dynamics of the membrane (action) potential V , Rubin and Wechselberger [24] consider that m and V evolve on the same fast time scale. So, the Hodgkin-Huxley model may be transformed into a singularly perturbed system with two time scales in which the slow variables are (n, h) and the fast variables are (m, V ).
Let's notice that the multiplicative parameter ε has been introduced artificially in Eq. (48c). This is due to the fact that it has been stated above that the time scale of variable m, i.e., y 1 is tenth times greater than the time scale of variables n and h, i.e. of variables x 1 and x 2 . Moreover, this parameter is identical to those use in Eq. (48d) since it has been also considered that m and V , i.e., y 1 and y 2 evolve on the same fast time scale.
According to the Geometric Singular Perturbation Theory, the zero-order approximation in ε of the slow manifold associated with the Hodgkin-Huxley model (48) is obtained by posing ε = 0 in Eqs. (48c & 48d). So, the slow manifold is given by: Then, the fast foliation is within the planes x 1 = constant and x 2 = constant.
The fold curve is defined as the location of the points where g 1 (x 1 , x 2 , y 1 , y 2 ) = 0, g 2 (x 1 , x 2 , y 1 , y 2 ) = 0 and det J (g1,g2) = 0. For the Hodgkin-Huxley model (48), the fold curve is thus given by Eqs. (49a & 49b) and by the determinant of the Jacobian matrix of the following fast foliation: The Jacobian matrix of the fast foliation (50) reads: where the ( ′ ) denotes the derivative with respect to y 2 . Then, taking into account Eqs. (49b), i.e., y 1 =ȳ 1 we have: So, the determinant of the Jacobian matrix of the fast foliation (50) is: Thus, the condition for the fold curve is det J (g1,g2) = 0, which gives: Therefore: By subtracting Eq. (49a) from Eq. (55) we obtain x 1 : Plugging this value of x 1 (56) into Eq. (55) provides: So, the fold curve is given by the set of parametric equations (56-57) in terms of y 2 .
The pseudo singular points are given by Eqs. (20) which reads for the Hodgkin-Huxley model (48): Let's notice that Eqs. (58c) and (58d) are identical. Moreover, the definition of τ 3 (46f) enables to simplify the above system (58). Thus, we have:Ī Moreover, Eqs. (59a) and (59c) indicate that the pseudo singular point belongs to the slow manifold and to the fold curve. So, let's replace in Eq. (59b) the variables x 1 and x 2 by the variables x 1f and x 2f given by Eq. (56) and Eq. (57) respectively which represent the parametric equations of fold curve.
Thus, it appears that Eq. (60) depends on the variable y 2 , on the functions gates dynamics time scales τ 1 (y 2 ) and τ 2 (y 2 ) and on the bifurcation parameterĪ. According to Rubin and Wechselberger [24], the function y 2 (Ī), solution of (60) is independent of time multiplicative constants k 1 and k 2 that one could set in factor of τ 1 (y 2 ) and τ 2 (y 2 ). So, following their works, let's plot the function y 2 (Ī) solution of (60) for various values of these time constants by posing successively in (60) k 1 = 1, 3, 4.75 and 7 and while fixing k 2 = 1. The result is presented in Fig. 2. Let's notice that this plot 8 is exactly the same as those presented by Rubin and Wechselberger [24] (Fig. 8-9) for a nondimensionalized threedimensional Hodgkin-Huxley singularly perturbed system which had been obtained after the following variable changes: V → −V andĪ → −I, then V → V + 65 and finally V → V /100.
We observe from Fig. 2 that the bifurcation parameter valueĪ C ≈ −4.8 is exactly identical (in absolute value) to those obtained by Rubin and Wechselberger [24]. Numerical resolution 9 of Eq. (60) provides a better approximation of the bifurcation parameter value: This value corresponds to a voltage y 2 = −3.18136 mV .
In Fig. 3, 4 & 5 canard solution of the four-dimensional Hodgkin-Huxley singularly perturbed system for the "canard value" ofĪ ≈ −4.1 has been plotted in the (x 1 , x 2 , y 2 ) phase-space and then in the (x 1 , y 2 ) phase plane. The green point represents the pseudo singular point. The trajectory curve, i.e., the canard solution has been plotted in red while the fold curve is in yellow. We observe on Fig. 3 that when the trajectory curve reaches the fold at the pseudo singular point it "jump" suddenly to the other part of the slow manifold before being reinjected towards the pseudo singular point.

DISCUSSION
In a previous paper entitled: "Canards Existence in Memristor's Circuits" (see Ginoux & Llibre [17]) 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. This method enabled to state a unique "generic" condition for the existence of "canard solutions" for such three and fourdimensional 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 [6] and then by Szmolyan and Wechselberger [25] and finally by Wechselberger [39] was considered as "generic" since it was exactly the same for singularly perturbed systems of dimension three and four with only one fast variable. In this work we have extended this new method to the case of four-dimensional singularly perturbed systems with two slow and two fast variables and we have stated that the condition for the existence of "canard solutions" in such systems is exactly identical to those proposed in our previous paper. This result confirms the genericity of the condition (σ 2 < 0) we have highlighted and provides a simple and efficient tool for testing the occurrence of "canard solutions" in any three or four-dimensional singularly perturbed systems with one or two fast variables. Applications of this method to the famous coupled FitzHugh-Nagumo equations and to the Hodgkin-Huxley model has enabled to show the existence of "canard solutions" in such systems. However, in this paper, only the case of pseudo singular points or folded singularities of saddle-type has been analyzed. Of course, the case of of pseudo singular points or folded singularities of node-type and focus-type could be also studied with the same method.

APPENDIX
Change of coordinates leading to the normal forms of four-dimensional singularly perturbed systems with two fast variables are given in the following section.