Lorenz-Krishnamurthy Model

The Slow Invariant Manifold of the

Keywords: Lorenz-Krishnamurthy model, Flow Curvature Method, Darboux invariance, Fenichel theory, slow invariant manifold

During this last decades, several attempts to construct slow invariant manifold of the Lorenz-Krishnamurthy five-mode model of slow-fast interactions in the atmosphere have been made by various authors. Unfortunately, as in the case of many two-time scales singularly perturbed dynamical systems the various asymptotic procedures involved for such a construction diverge. So, it seems that till now only the first-order and third-order approximations of this slow manifold have been analytically obtained. While using the Flow Curvature Method we show in this work that one can provide the eighteenth-order approximation of the slow manifold of the generalized Lorenz-Krishnamurthy model and the thirteenth-order approximation of the "conservative" Lorenz-Krishnamurthy model. The invariance of each slow manifold is then established according to Darboux invariance theorem.

Introduction

The classical geometric theory developed originally by Andronov [START_REF] Andronov | Theory of Oscillators, I[END_REF], Tikhonov [START_REF] Tikhonov | On the dependence of solutions of differential equations on a small parameter[END_REF] and Levinson [START_REF] Levinson | A second-order differential equation with singular solutions[END_REF] stated that singularly perturbed systems possess invari-Preprint submitted to Elsevier Science 13 August 2018

ant manifolds on which trajectories evolve slowly and toward which nearby orbits contract exponentially in time (either forward and backward) in the normal directions. These manifolds have been called asymptotically stable (or unstable) slow manifolds. Then, Fenichel [START_REF] Fenichel | Persistence and Smoothness of Invariant Manifolds for Flows[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] theory 1 for the persistence of normally hyperbolic invariant manifolds enabled to establish the local invariance of slow manifolds that possess both expanding and contracting directions and which were labeled slow invariant manifolds. During the last century, various methods have been developed in order to determine the slow invariant manifold analytical equation associated with singularly perturbed systems. The seminal works of Wasow [START_REF] Wasow | Asymptotic Expansions for Ordinary Differential Equations[END_REF], Cole [START_REF] Cole | Perturbation Methods in Applied Mathematics[END_REF], O'Malley [START_REF] O'malley | Introduction to Singular Perturbations[END_REF][START_REF] O'malley | Singular Perturbations Methods for Ordinary Differential Equations[END_REF] and Fenichel [START_REF] Fenichel | Persistence and Smoothness of Invariant Manifolds for Flows[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] to name but a few, gave rise to the so-called Geometric Singular Perturbation

Method. According to this theory, existence as well as local invariance of the slow manifold of singularly perturbed systems have been stated.

Then, the determination of the slow manifold analytical equation turned into a regular perturbation problem [30, p. 112] in which one generally expected the asymptotic validity of such expansion to breakdown [START_REF] O'malley | Singular Perturbations Methods for Ordinary Differential Equations[END_REF].

systems of ordinary differential equations with two time scales, called Flow Curvature Method [START_REF] Ginoux | Slow manifold of a neuronal bursting model[END_REF][START_REF] Ginoux | Differential Geometry and Mechanics Applications to Chaotic Dynamical Systems[END_REF][START_REF] Ginoux | Slow Invariant Manifolds as Curvature of the Flow of Dynamical Systems[END_REF][START_REF] Ginoux | Differential geometry applied to dynamical systems[END_REF][START_REF] Ginoux | The flow curvature method applied to canard explosion[END_REF][START_REF] Ginoux | Canards from Chua's circuit[END_REF] has been developed. It consists in considering the trajectory curves integral of such systems as curves in Euclidean n-space.

Based on the use of local metrics properties of curvatures inherent to Differential Geometry, this method which does not require the use of asymptotic expansions, states that the location of the points where the local curvature of the trajectory curves of such systems is null defines a (n -1)-dimensional manifold associated with this system and called flow curvature manifold. The invariance of this manifold is then stated according to a theorem introduced by Gaston Darboux [START_REF] Darboux | Mmoire sur les quations diffrentielles algbriques du premier ordre et du premier degr[END_REF] in 1878.

The laws governing the behavior of the atmosphere permit the simultaneous presence of oscillation modes such as quasi-geostrophic modes and inertialgravity modes. The former which have periods of few days are generally referred as Rossby Waves while the latter whose periods are of few hours are called Gravity Waves. In 1980, a nine-dimensional primitive equation (PE) model of the atmosphere enabling to superpose Rossby and Gravity Waves was originally proposed by the late Edward Norton Lorenz [START_REF] Lorenz | Attractor sets and quasi-geostrophic equilibrium[END_REF]. Few years later, Lorenz [START_REF] Lorenz | On the existence of a slow manifold[END_REF] simplified the nine-dimensional model to a five-dimensional model by truncating to just five modes: three Rossby Waves coupled to two Gravity Waves. This five-dimensional model can be considered as a two-time scales singularly perturbed dynamical system with three slow variables (Rossby Waves) and two fast variables (Gravity Waves) in accordance with physical observations of atmospheric behavior. In numerical weather prediction a problem arises because raw fields data can not be used as initial conditions for such model, since even when the initial wind and pressure fields are both fairly re-alistic, Gravity Waves will occur, if the fields are not in "proper balance".

According to Camassa [3, p. 357]: "Small errors in the "proper balance" between these two time scales lead to abnormal evolution of gravity waves, which in turn causes appreciable deviation of weather forecasts from actual observation on the time scale of gravity waves.". To solve this initialization problem, existence of a slow manifold2 , consisting of trajectory curves (orbits) for which Gravity Waves motion is absent, was first postulated for such model. Then, an iteration scheme was developed to find from the state (point in phase space)

specified by field data a corresponding initial state on this slow manifold, so that weather forecasts with these initial states can be accurate on the same time scale as Rossby Waves. In their paper, Lorenz and Krishnamruthy [START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF] identified the variables representing Gravity Wave activity as the ones which can exhibit fast oscillations, and defined the slow manifold as an invariant manifold in the five dimensional phase space for which fast oscillations never develop. However, in a subsequent paper, Lorenz and Krishnamruthy [START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF] identified a trajectory curve (orbit) which by construction has to lie on the slow manifold. They followed its evolution numerically to show that sooner or later fast oscillations developed, thereby implying, as pointed out by the title of their article, that a slow manifold according to the definition does not exist for this model. Such a result gave rise to a series of articles published from 1991 to 1996 by S.J. Jacobs [START_REF] Jacobs | Existence of a slow manifold in a model system of equations[END_REF], J.P. Boyd [START_REF] Boyd | The slow manifold of five-mode model[END_REF], A.C. Fowler and G. Kember [START_REF] Fowler | The LorenzKrishnamurthy slow manifold[END_REF],

R. Camassa and Siu-Kei Tin [START_REF] Camassa | The global geometry of the slow manifold in the Lorenz Krishnamurthy model[END_REF] in which their authors proved the existence of a slow manifold for the Lorenz-Krishnamurthy model (LK). These latter proposed a generalized LK model. More recently, M. Phani Sudheer, Ravi S.

Nanjundiah A.S. Vasudeva Murthy [START_REF] Sudheer | Revisiting the slow manifold of the Lorenz-Krishnamurthy quintet[END_REF] "revisited" the same model and pro-vided an approximation of its slow manifold while using Girimaji's technique of local reduction. Then, J. Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] studied a "conservative form" of the Lorenz-Krishnamurthy model and developed distinct methods to derive the leading-order asymptotics of the late coefficients in the power-series expansion used for the construction of its slow manifold. But, to our knowledge, only the first-order and third-order approximations of this slow manifold have been analytically obtained. So, the aim of this paper is to show that while using the Following the approach of C.K.R.T. Jones [START_REF] Jones | Geometric Singular Pertubation Theory[END_REF] and Kaper [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] some fundamental concepts and definitions for systems of ordinary differential equations with two time scales, i.e., for singularly perturbed dynamical systems are briefly recalled.

In the following we consider a dynamical systems theory for systems of differential equations of the form:

               x = f ( x, z, ε) z = ε g ( x, z, ε) (1)
where In the case when ε 1, i.e., is a small positive number, the variable x is called fast variable, and z is called slow variable. Using Landau's notation:

x ∈ R m , z ∈ R p , ε ∈ R + ,
O(ε l ) represents a real polynomial in ε of l degree, with l ∈ Z, it is used to consider that generally x evolves at an O (1) rate; while z evolves at an O (ε) slow rate.

Reformulating the system (1) in terms of the rescaled variable τ = εt, we obtain:

               ε ˙ x = f ( x, z, ε) ˙ z = g ( x, z, ε) (2)
The dot (•) represents the derivative with respect to the new independent variable τ .

The independent variables t 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 ε = 0, and they are labeled singular perturbation problems when ε 1, i.e., is a small positive parameter. The label "singular" stems in part from the discontinuous limiting behavior in the system (1) as ε → 0.

In such case, the system (1) reduces to an m-dimensional system called reduced fast system, with the variable z as a constant parameter:

               x = f ( x, z, 0) z = 0 (3)
System (2) leads to the following differential-algebraic system called reduced slow system which dimension decreases from m + p to p:

               0 = f ( x, z, 0) ˙ z = g ( x, z, 0) (4)
By exploiting the decomposition into fast and slow reduced systems (3) and ( 4), the geometric approach reduced the full singularly perturbed system to separate lower-dimensional regular perturbation problems in the fast and slow regimes, respectively.

Fenichel geometric theory

Fenichel geometric theory for general systems (1), i.e., a theorem providing conditions under which normally hyperbolic invariant manifolds in system [START_REF] Andronov | Theory of Oscillators, I[END_REF] persist when the perturbation is turned on, i.e., when 0 < ε 1 is briefly recalled in this subsection. This theorem concerns only compact manifolds with boundary.

Normally hyperbolic manifolds

Let's make the following assumptions about system (1): (H 3 ) M 0 is normally hyperbolic relative to (3) and in particular it is required for all points p ∈ M 0 , that there are k (resp. l) eigenvalues of D y f ( p, 0) with positive (resp. negative) real parts bounded away from zero, where k + l = m.

(H 1 )

Fenichel persistence theory for singularly perturbed systems

For compact manifolds with boundary, Fenichel's persistence theory states that, provided the hypotheses (H 1 ) -(H 3 ) are satisfied, the system (1) has a slow (or center) manifold, and this slow manifold has fast stable and unstable manifolds.

Theorem for compact manifolds with boundary:

Let system (1) satisfy the conditions (H Proof. For proof of this theorem see Fenichel [START_REF] Fenichel | Persistence and Smoothness of Invariant Manifolds for Flows[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF].

The label slow manifold is attached to M ε because the magnitude of the vector field restricted to M ε is O (ε), in terms of the fast independent variable t.

So persistent manifolds are labeled slow manifolds, and the proof of their persistence is carried out by demonstrating that the local stable and unstable manifolds of M 0 also persist as locally invariant manifolds in the perturbed system, i.e., that the local hyperbolic structure persists, and then the slow manifold is immediately at hand as a locally invariant manifold in the transverse intersection of these persistent local stable and unstable manifolds.

Geometric singular perturbation theory

Earliest geometric approaches to singularly perturbed systems have been developed by Cole [START_REF] Cole | Perturbation Methods in Applied Mathematics[END_REF], O'Malley [START_REF] O'malley | Introduction to Singular Perturbations[END_REF][START_REF] O'malley | Singular Perturbations Methods for Ordinary Differential Equations[END_REF], Fenichel [START_REF] Fenichel | Persistence and Smoothness of Invariant Manifolds for Flows[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions[END_REF][START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] for the determination of the slow manifold equation.

Generally, Fenichel theory enables to turn the problem for explicitly finding functions y = Y ( x, ε) whose graphs are locally invariant slow manifolds M ε of system (1) into regular perturbation problem [30, p. 112]. Invariance of the manifold M ε implies that Y ( x, ε) satisfies:

D x Y ( x, ε) f x, Y ( x, ε) , ε = ε g x, Y ( x, ε) , ε (5)
According to Guckenheimer et al. [18, p. 131], this (partial) differential equation for Y ( x, ε) cannot be solved exactly. So, its solution can be approximated arbitrarily closely as a Taylor series at ( x, ε) = 0, 0 .

Then, the following perturbation expansion is plugged:

Y ( x, ε) = Y 0 ( x) + ε Y 1 ( x) + O ε 2 (6)
into [START_REF] Cole | Perturbation Methods in Applied Mathematics[END_REF] to solve order by order for Y ( x, ε). The Taylor series expansion [START_REF] O'malley | Singular Perturbations Methods for Ordinary Differential Equations[END_REF] for f x, Y ( x, ε) , ε and g x, Y ( x, ε) , ε up to terms of order two in ε reads:

f x, Y ( x, ε) , ε = f x, Y 0 ( x) , 0 + ε   D y f x, Y 0 ( x) , 0 Y 1 ( x) + ∂ f ∂ε x, Y 0 ( x) , 0   g x, Y ( x, ε) , ε = g x, Y 0 ( x) , 0 + ε D y g x, Y 0 ( x) , 0 Y 1 ( x) + ∂ g ∂ε x, Y 0 ( x) , 0
• At order ε 0 , Eq. ( 5) gives:

D x Y 0 ( x) f ( x, Y 0 ( x) , 0) = 0 (7)
which defines Y 0 ( x) due to the invertibility of D y f and the Implicit Function Theorem.

• The next order ε 1 provides:

D x Y 0 ( x)   D y f ( x, Y 0 ( x) , 0) Y 1 ( x) + ∂ f ∂ε   = g( x, Y 0 ( x) , 0) (8)
which yields Y 1 ( x) and so forth.

So, regular perturbation theory makes it possible to build an approximation of locally invariant slow manifolds M ε . Thus, in the framework of the Geometric Singular Perturbation Method, three conditions are needed to characterize the slow manifold associated with singularly perturbed system: existence, local invariance and determination. Existence and local invariance of the slow manifold are stated according to Fenichel theorem for compact manifolds with boundary while asymptotic expansions provide its equation up to the order of the expansion.

Flow Curvature Method

In this section, one of the main results of the Flow Curvature Method and based on the use of local properties of curvatures inherent to Differential Geometry is briefly presented (for more details see [START_REF] Ginoux | Slow manifold of a neuronal bursting model[END_REF][START_REF] Ginoux | The flow curvature method applied to canard explosion[END_REF]). According to this method, the highest curvature of the flow, i.e. the (n-1) th curvature of trajectory curve integral of n-dimensional singularly perturbed dynamical systems (1) defines a (n -1)-dimensional manifold associated with this system and called flow curvature manifold. We have the following result:

Slow manifold equation

Proposition 1 The location of the points where the (n -1) th curvature of the flow, i.e. the curvature of the trajectory curve X, integral of any n-dimensional singularly perturbed dynamical systems (1) vanishes, provides a k-order approximation in ε of its slow manifold M ε the equation of which reads

φ( X, ε) = ˙ X • ( ¨ X ∧ ... X ∧ . . . ∧ (n) X ) = det( ˙ X, ¨ X, ... X, . . . , (n)
X ) = 0 (9) where (n)
X represents the time derivatives up to order n of X = ( x, z) t .

Note.

k-order approximation depends on the number of ε contained in the vector field. We will see in Sec. 4 that for the Lorenz-Krishnamurthy model k = 18.

Proof.

While the slow invariant manifold analytical equation [START_REF] Darboux | Mmoire sur les quations diffrentielles algbriques du premier ordre et du premier degr[END_REF] given by the Geometric Singular Perturbation Method is an explicit equation, the slow invariant manifold analytical equation ( 9) obtained according to the Flow Curvature Method is an implicit equation. So, in order to compare the latter with the former it is necessary to plug the following perturbation expansion: 9). Thus, solving order by order for Y ( x, ε) will transform ( 9) into an explicit analytical equation enabling the comparison with [START_REF] Darboux | Mmoire sur les quations diffrentielles algbriques du premier ordre et du premier degr[END_REF]. The Taylor series expansion for φ( X, ε) = φ( x, Y ( x, ε) , ε) up to terms of order one in ε reads:

Y ( x, ε) = Y 0 ( x) + ε Y 1 ( x) + O (ε 2 ) into (
φ( X, ε) = φ( x, Y 0 ( x) , 0) + εD y φ( x, Y 0 ( x) , 0) Y 1 ( x) + ε ∂φ ∂ε ( x, Y 0 ( x) , 0) (10)
• At order ε 0 , Eq. ( 10) gives:

φ x, Y 0 ( x) , 0 = 0 (11)
which defines Y 0 ( x) due to the invertibility of D y φ and application of the Implicit Function Theorem.

• The next order ε 1 , provides:

D y φ( x, Y 0 ( x) , 0) Y 1 ( x) + ∂φ ∂ε ( x, Y 0 ( x) , 0) = 0 ( 12
)
which yields Y 1 ( x) and so forth.

In order to prove that this equation is completely identical to Eq. ( 8), let's rewrite it as follows:

Y 1 ( x) = -D y φ( x, Y 0 ( x) , 0) -1 ∂φ ∂ε ( x, Y 0 ( x) , 0)
By application of the chain rule, i.e., the derivative of φ( x, Y 0 ( x) , 0) with respect to the variable y and then with respect to ε, it can be stated that:

Y 1 ( x) = -(D x f )(D y f ) -1 (D y f ) g( x, Y 0 ( x), 0) -D y f -1 D ε f ( x, Y 0 ( x) , 0)
But, according to the Implicit Function Theorem we have:

(D x f ) = -(D y f )(D x y) = -(D y f )(D x Y 0 ( x))
Then, by replacing into the previous equation we find:

Y 1 ( x) = (D y f )(D x Y 0 ( x))(D y f ) -1 (D y f ) g( x, Y 0 ( x), 0)-D y f -1 D ε f ( x, Y 0 ( x) , 0)
After simplifications, we have:

Y 1 ( x) = D x Y 0 ( x) D y f ( x, Y 0 ( x) , 0) -1 g( x, Y 0 ( x) , 0)-D y f ( x, Y 0 ( x) , 0) -1 ∂ f ∂ε 14
Finally, Eq. ( 12) may be written as:

D x Y 0 ( x)   D y f ( x, Y 0 ( x) , 0) Y 1 ( x) + ∂ f ∂ε   = g( x, Y 0 ( x) , 0)
Thus, identity between the "slow manifold " equation given by the Geometric Singular Perturbation Method and by the Flow Curvature Method is proved up to first order term in ε.

Note. Let's notice that the slow invariant manifold equation ( 9) associated with n-dimensional singularly perturbed systems defined by the Flow Curvature Method is a tensor of order n. As a consequence, it can only provide an approximation of n-order in ε of the slow invariant manifold equation ( 6).

Nevertheless, it is easy to show that the Lie derivative of the "slow manifold " equation ( 9) obtained by the Flow Curvature Method can be written as:

L X φ( X, ε) = ˙ X • ( ¨ X ∧ ... X ∧ . . . ∧ (n+1) X ) = det( ˙ X, ¨ X, ... X, . . . , (n+1)
X ) = 0 [START_REF] Ginoux | Differential Geometry and Mechanics Applications to Chaotic Dynamical Systems[END_REF] where

(n+1)
X represents the time derivatives up to order (n + 1) of X = ( x, y) t . So, Eq. ( 13) defines a tensor of order n+1 which provides an approximation of (n + 1)-order in ε of the slow invariant manifold equation ( 6). Thus, by taking the successive Lie derivatives of the "slow manifold " equation ( 9) we improve the order of the approximation up to an order corresponding to that of the Lie derivative. As an example, according to Prop. 1, the "slow manifold " equation of a two-dimensional singularly perturbed dynamical system reads:

φ( X, ε) = det( ˙ X, ¨ X) = 0
where X = (x, y) t . This second-order tensor only provides a first order approximation in ε of the slow invariant manifold equation ( 6). While its Lie derivative

L X φ( X, ε) = det( ˙ X, ... X) = 0
which is third-order tensor gives a second-order approximation in ε. Thus, by applying a mathematical induction, the proof above can be extended to high-order approximations in ε.

Invariance of the slow manifold

The local invariance of the slow manifold analytical equation defined by Flow Curvature Method may be stated while using the Tangent Linear System Approximation (T.L.S.A.) associated with Darboux Invariance Theorem. Tangent Linear System Approximation has been introduced by Rossetto et al. [33] in order to compute the slow manifold equation of singularly perturbed systems. This approximation consists in replacing, in the vicinity of the singular approximation, the singularly perturbed system by the corresponding tangent linear system. Tangent Linear System Approximation may thus be viewed as a local invariance condition of the slow manifold.

Tangent linear system approximation (T.L.S.A.)

The tangent linear system approximation (T.L.S.A.) states that, in the vicinity of the singular approximation associated with singularly perturbed systems (1), the functional jacobian matrix of such systems is locally stationary, i.e., dJ dt = 0 (14)

Lie Derivative -Darboux Invariance Theorem

Let φ a C 1 function defined in a compact E included in R and X (t) the integral of the dynamical system defined by [START_REF] Andronov | Theory of Oscillators, I[END_REF]. The Lie derivative is defined as follows:

L X φ = ˙ X • - → ∇φ = n i=1 ∂φ ∂x i ẋi = dφ dt (15)
Darboux Invariance Theorem:

An invariant manifold is defined by φ( X, ε) = 0 where φ is a C 1 in an open set U and such there exists a C 1 function denoted k( X) and called cofactor which satisfies

L X φ( X, ε) = k( X)φ( X, ε) for all X ∈ U (16)
Proof.

The proof of this theorem is in Darboux [START_REF] Darboux | Mmoire sur les quations diffrentielles algbriques du premier ordre et du premier degr[END_REF]. However, let's prove that both Darboux and Fenichel's invariance are exactly identical.

According to Fenichel's persistence theorem a slow invariant manifold M ε may be written as an explicit function: y = Y ( x, ε) the invariance of which implies that Y ( x, ε) satisfies:

D x Y ( x, ε) f ( x, Y ( x, ε) , ε) = ε g( x, Y ( x, ε) , ε) (17)
Let's write the slow manifold M ε as an implicit function by posing:

φ( x, y, ε) = y -Y ( x, ε) (18)
According to Darboux Invariance Theorem M ε is invariant if and only if:

L V φ( x, y, ε) = k( x, y, ε)φ( x, y, ε) (19)
Plugging Eq. ( 18) into the Lie derivative ( 19) leads to:

L V φ( x, y, ε) = ˙ y -D x Y ( x, ε) ˙ x = k( x, y, ε)φ( x, y, ε)
which may be written according to Eq. ( 1):

L V φ( x, y, ε) = ε g( x, y, ε) -D x Y ( x, ε) f ( x, y, ε) = k( x, y , ε)φ( x, y, ε)
Evaluating this Lie derivative in the location of the points where φ( x, y, ε) = 0, i.e. y = Y ( x, ε) leads to:

L V φ( x, Y ( x, ε) , ε) = ε g( x, Y ( x, ε) , ε) -D x Y ( x, ε) f ( x, Y ( x, ε) , ε) = 0
which is exactly identical to Eq. ( 5) proposed by Fenichel.

L V φ( X) = ˙ X • ( ¨ X ∧ ... X ∧ . . . ∧ J (n) X ) (24)
The right hand side of this Eq. ( 24) can be written:

J ˙ X•( ¨ X ∧ ... X ∧ . . . ∧ (n) X )+ ˙ X•(J ¨ X ∧ ... X ∧ . . . ∧ (n) X )+. . .+ ˙ X•( ¨ X ∧ ... X ∧ . . . ∧ J (n) X )
According to Eq. ( 23) all terms are null except the last one. So, by taking into account identity (42) established in Appendix we find:

L V φ( X) = T r (J) ˙ X • ( ¨ X ∧ ... X ∧ . . . ∧ (n) X ) = T r (J) φ( X) = k( X)φ( X)
where k( X) = T r (J) represents the trace of the functional jacobian matrix.

So, according to Prop. 2 invariance of the slow manifold analytical equation of any n-dimensional singularly perturbed dynamical system is established provided that the functional jacobian matrix is locally stationary.

The Lorenz-Krishnamurthy Slow Invariant Manifold

Let's consider the model introduced by the late E.N. Lorenz [START_REF] Lorenz | On the existence of a slow manifold[END_REF] and usually referred to as the Lorenz five-mode model [START_REF] Boyd | The slow manifold of five-mode model[END_REF] or as the Lorenz-Krishnamurthy model [START_REF] Lorenz | On the existence of a slow manifold[END_REF][START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF][START_REF] Lorenz | The slow manifold What is it?[END_REF]. This model, obtained by truncation of the rotating shallowwater equations, governs the dynamics of a triad of vortical modes, with amplitudes (u, v, w), coupled to a gravity mode described by (x, y). In 1996, Camassa et al. [START_REF] Camassa | The global geometry of the slow manifold in the Lorenz Krishnamurthy model[END_REF] proposed a generalized LK model presented in the next section.

                                                   δ ẋ = -y -κx δ ẏ = x + uv -κy u = -vw + δ vy -δu v = uw -δ uy -δv + F ẇ = -(uv + δw) (27)
In the first article in which he introduced his five-mode model, Lorenz [23, p. 1548] wrote: "In the trivial case where the Rossby and gravity waves are completely uncoupled, i.e., where the system of equations degenerates into two systems, one governing Rossby waves and one governing gravity waves, the slow manifold obviously exists and is obtained simply by equating all of the gravity-waves variables to zero."

Following this idea Camassa et al. [4, p. 3263] gave the zero-order approximation in δ of the slow manifold associated with the system (27):

x = - uv 1 + κ 2 + O(δ) y = κ uv 1 + κ 2 + O(δ) (28)
This slow manifold parametrized as a graph over the (u, v, w) space is also called singular approximation (see Sec. 2 for definition) since it is obtained by posing δ = 0 in the two first equations of system [START_REF] Lorenz | On the existence of a slow manifold[END_REF].

After the publication of the article entitled "On the Nonexistence of a Slow Manifold" in which Lorenz and Krishnamurthy [START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF] concluded that the slow manifold of such model "does not exist", Jacobs [START_REF] Jacobs | Existence of a slow manifold in a model system of equations[END_REF], Boyd [START_REF] Boyd | The slow manifold of five-mode model[END_REF], Fowler and Kember [START_REF] Jacobs | Existence of a slow manifold in a model system of equations[END_REF] and then Camassa and Tin [START_REF] Camassa | The global geometry of the slow manifold in the Lorenz Krishnamurthy model[END_REF] proved the existence of a slow manifold in this model and gave approximations of its equation at first orders while using the "singular perturbation scheme known as the method of multiple scales. 4 ". Although these authors stated that a slow manifold can be constructed via formal series, such a long and tedious asymptotic procedure of systematic identification order-by-order is expected to diverge as previously recalled. Recently, Sudheer et al. [START_REF] Sudheer | Revisiting the slow manifold of the Lorenz-Krishnamurthy quintet[END_REF] and Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] have proposed alternative techniques for the construction of the slow manifold of such model.

We will show now that one can obtain high-orders approximation of this slow manifold while using the Flow Curvature Method.

Thus, according to Prop. 1 the slow manifold equation ( 9) associated with the generalized LK model (27) reads5 :

φ( X, δ) = det( ˙ X, ¨ X, ... X, .... X , ..... X ) = 0 (29)
where X = ( x, z) t . Then, it can be verified that the time derivative of the functional jacobian matrix of system (26) evaluated when δ → 0 is a zero matrix. So, from Darboux Invariance Theorem we can conclude that in the δ-vicinity of the singular approximation the slow manifold is invariant.

The implicit equation ( 29) is a polynomial of degree 10 for u, v and w, of degree 5 for x and 11 for y and represents the eighteenth-order approximation in δ of the slow manifold of the generalized LK model [START_REF] Lorenz | On the existence of a slow manifold[END_REF].

By posing δ = 0 in the above Eq. ( 29) we find that:

φ( X, 0) = v(u 2 -w 2 ) (x + uv 1 + κ 2 ) 2 + (y -κ uv 1 + κ 2 ) 2 = 0 ( 30
)
This equation is made of a product of invariant manifolds as it is easy to verify according to Darboux Invariance Theorem. Let's compute the Lie derivative of the first and second term, when δ → 0 we have:

L X (v) = 0 L X (u 2 -w 2 ) = 0 (31)
Let's notice that the third term of Eq. ( 30) is nothing else but the zero-order approximation in δ (singular approximation) of the slow manifold (see Eq.

(

)) given by Camassa et al. [4, p. 3263] which is also invariant when δ → 0.

But, according to Leith [24,p. 960], the decomposition into fast and slow modes enables to define a three-dimensional submanifold of the state space parametrized by (u, v, w) and that he called slow manifold. So, let's pose x → 0 and y → 0 in the above Eq. ( 29) we find that: As previously, it can be verified that the time derivative of the functional jacobian matrix of the fast system (26) (from which the slow system (34) has been deduced) is a zero matrix when ε → 0. So, from Darboux Invariance

Theorem we can conclude that in the ε-vicinity of the singular approximation the slow manifold is invariant.

The implicit equation ( 37) is a polynomial of degree 9 for u, v and w, of degree 5 for x and 11 for y and represents the thirteenth-order approximation in ε of the slow manifold of the conservative LK model [START_REF] Tikhonov | On the dependence of solutions of differential equations on a small parameter[END_REF].

By posing ε = 0 in the above Eq. ( 37) we find that:

φ( X, 0) = (u 2 -w 2 )(v 2 + w 2 )((x + buv) 2 + y 2 ) = 0 (38)

This slow manifold is made of a product of invariant manifolds as it is easy to verify according to Darboux Invariance Theorem. Let's compute the Lie derivative of the first and second term when ε → 0 we have: model (34) has been plotted in Fig. 1 in the (u, v, w) phase-space. Numerical integration of this model with a set of initial conditions (x 0 , y 0 , u 0 , v 0 , w 0 ) = (2, 2, -2, 1.97, 2) taken on this slow manifold (in blue on Fig. 1) enables to highlight that the trajectory curves (in red on Fig. 1) "visit" every part of this hypersurface and stay in its ε-vicinity. The fixed point located at the origin has been plotted in green in the center of this figure.

  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, z and ε in U ×I, where U is an open subset of R m ×R p and I is an open interval containing ε = 0.

  The functions f and g are C ∞ in U × I, where U is an open subset of R m × R p and I is an open interval containing ε = 0. (H 2 ) There exists a set M 0 that is contained in ( x, y) : f ( x, y, 0) = 0 such that M 0 is a compact manifold with boundary and M 0 is given by the graph of a C ∞ function y = Y 0 ( x) for x ∈ D, where D ⊆ R p is a compact, simply connected domain and the boundary of D is an (p -1) dimensional C ∞ submanifold. Finally, the set D is overflowing invariant with respect to (2) when ε = 0.

  [START_REF] Andronov | Theory of Oscillators, I[END_REF] ) -(H 3 ). If ε > 0 is sufficiently small, then there exists a function Y ( x, ε) defined on D such that the manifold M ε = ( x, y) : y = Y ( x, ε) is locally invariant under[START_REF] Andronov | Theory of Oscillators, I[END_REF]. Moreover, Y ( x, ε) is C r for any r < +∞, and M ε is C r O (ε) close to M 0 . In addition, there exist perturbed local stable and unstable manifolds of M ε . They are unions of invariant families of stable and unstable fibers of dimensions l and k, respectively, and they are C r O (ε) close for all r < +∞ , to their counterparts.

φ 2 ( 2 + 2 )

 222 (u, v, w, δ) = u 2 w 2 (v 2 δ 2 (1 + 2 ) -(δ 2 -κ) 2 ) -w 2 (u 2 + v 2 w 2 δ 2 ) + u 4 (1 + (δ 2 -κ) 2 ) + uvwδ(δ 2 -κ)(2w2 -u using the Flow Curvature Method, i.e. according to Prop. 1 the slow manifold equation (9) associated with conservative LK model (34) reads 6 :

Fig. 1 .

 1 Fig. 1. The conservative LK model slow invariant manifold in (u, v, w)-space

This theory was independently established in Hirsch, et al.[START_REF] Hirsch | Invariant Manifolds[END_REF]

This concept has been introduced by C. E. Leith[START_REF] Leith | Nonlinear Normal Mode Initialization and Quasi-Geostrophic Theory[END_REF] in 1980.

In certain applications these functions will be supposed to be C r , r 1

Boyd [2, p. 1058]

See http://ginoux.univ-tln.fr for complete equation.

See http://ginoux.univ-tln.fr for complete equation

Acknowledgements

Authors would like to thank Prof. Jaume Llibre for his mathematical remarks and comments.

Now, let's prove the invariance of the flow curvature manifold [START_REF] Fenichel | Asymptotic stability with rate conditions II[END_REF], i.e., the invariance of the "slow manifold equation" defined by Prop. 1.

Proposition 2

The flow curvature manifold defined by φ( X) = 0 where φ is a C 1 in an open set U is invariant with respect to the flow of (1) if there exists a C 1 function denoted k( X) and called cofactor which satisfies:

for all X ∈ U and where

Proof. Lie derivative of the flow curvature manifold (9) reads:

From the identity ¨ X = J ˙ X where J is the functional jacobian matrix associated with any n-dimensional singularly perturbed system (1) we find that:

where J n represents the n th power of J, e.g., ¨

Then, it follows that:

X in Eq. ( 14) with Eq. ( 16) we have:

Starting from this model we will provide the eighteenth-order approximation of its slow manifold while using the Flow Curvature Method, the invariance of which will be established according to Darboux theroem. Moreover, by posing δ = 0 in our slow manifold analytical equation we will find again the first-order approximation of the slow manifold given by Camassa et al. [4, p. 3263].

The generalized LK model

According to Camassa et al. [START_REF] Camassa | The global geometry of the slow manifold in the Lorenz Krishnamurthy model[END_REF], the Lorenz-Krishnamurthy model [START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF] can be written as:

where 0 is the coupling parameter between the Rossby Wave (u, v, w) and Gravity Wave (x, y), α is a parameter introduced to model the dissipation and controls the damping of Rossby mode, while κ controls the damping of the gravity mode. F > 0 represents the forcing parameter which is assumed to be much smaller than unity. By posing α = a, = b and κ = α in Eqs. [START_REF] Llibre | On the analytic integrability of the 5dimensional Lorenz system for the gravity-wave activity[END_REF] one finds again the Lorenz-Krishnamurthy model [START_REF] Lorenz | On the nonexistence of a slow manifold[END_REF]. Then, by making the following variables changing:

x → δ 2 x, y → δ 2 y, u → δu, v → δv, w → δw and while posing for convenience α = δ 2 , Camassa et al. [START_REF] Camassa | The global geometry of the slow manifold in the Lorenz Krishnamurthy model[END_REF] obtained the following system:

Note.

Reformulating the fast system (26) in terms of the rescaled variable τ = δt we obtain a two-time scales singularly perturbed dynamical system of the same form as the slow system (2) for which variables x = (x, y) t and z = (u, v, w) t are respectively fast and slow.

The conservative LK model

Thus, by using the same variables changing as previously and by posing δ = ε and b = in system [START_REF] Lorenz | On the existence of a slow manifold[END_REF], Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] obtained the following two-time scales singularly perturbed dynamical system:

where parameters b and ε control the strength of the coupling and the gravitywave frequency, x and y are fast modes while u, v and w are the slow modes.

At the beginning of his paper Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] gives the zero-order approximation in ε (singular approximation) of the slow manifold associated with the conservative LK model (34) by posing ε = 0:

Thus, as previously noticed by Lorenz [START_REF] Lorenz | On the existence of a slow manifold[END_REF] and recalled by Camassa [START_REF] Camassa | On the geometry of an atmospheric slow manifold[END_REF] and

Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] this model has an invariant manifold the equation of which is:

Let's notice that the third term of Eq. ( 38) is nothing else but the zero-order approximation in ε (singular approximation) of the slow manifold (see Eq. ( 35)) given by Vanneste [START_REF] Vanneste | Asymptotics of a slow manifold[END_REF] which is also invariant when ε → 0.

As previously, the decomposition into fast and slow modes enables to define a three-dimensional submanifold of the state space parametrized by (u, v, w)

and that he called slow manifold. So, let's pose x → 0 and y → 0 in the above Eq. ( 37) we find that:

In addition to the quadratic invariant manifolds (36-39) highlighted above we find another manifold. Let's compute its Lie derivative when ε → 0 we obtain:

Thus, we deduce that this manifold is locally invariant, i.e. is invariant in the vicinity of the manifold defined by u 2 -w 2 = 0.

APPENDIX

The identity involved in the proof of the invariance of the slow manifold (Sec.

4.2.2) is stated in this appendix.

Proof. The proof is based on inner product properties.

To the functional jacobian matrix J is associated an eigenbasis:

Let suppose that there exists a transformation 7 such that:

to each vector a i corresponds the eigenvector Y λ i with i = 1, . . . , n.

Each inner product of the left hand side Eq. ( 42) may be transformed into Making the sum of these factors the proof is stated.

7 By considering that each vector a i may be spanned on the eigenbasis, calculus is longer and tedious but leads to the same result.


Jean-Marc Ginoux
email: ginoux@univ-tln.fr

Jaume Llibre
email: jllibre@mat.uab.es

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

Keywords: Geometric singular perturbation theory, singularly perturbed dynamical systems, canard solutions

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 wellknown 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 [START_REF] Benoît | Les canards de R 3[END_REF], Benoît [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then Benoît [START_REF] Benoît | Canards de R 3[END_REF] 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 [START_REF] Szmolyan | Canards in R 3[END_REF] extended "Geometric Singular Perturbation Theory 2 " to canards problems in R 3 and provided a "standard version" of Benoît's theorem [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF]. Very recently, Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] 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 [START_REF] Szmolyan | Canards in R 3[END_REF] and Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] 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 and cj 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 [START_REF] Szmolyan | Canards in R 3[END_REF] nor Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] 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 [START_REF] Ginoux | Canards in Memristor's Circuits, Qualitative Theory of Dynamical Systems[END_REF]) 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 [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and finally by Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] 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 [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] 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.

DEFINITIONS 2.1. Singularly perturbed systems

According to Tikhonov [START_REF] Tikhonov | On the dependence of solutions of differential equations on a small parameter[END_REF], Jones [START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF] and Kaper [START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] singularly perturbed systems are defined as:

⃗ x ′ = ε ⃗ f (⃗ x, ⃗ y, ε) , ⃗ y ′ = ⃗ g (⃗ x, ⃗ y, ε) . ( 1
)
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 ∞ functions3 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

⃗ x = ⃗ f (⃗ x, ⃗ y, ε) , ε ⃗ y = ⃗ g (⃗ x, ⃗ y, ε) .
(

) 2
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 is

⃗ x = ⃗ f (⃗ x, ⃗ y, ε) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, ε) .
(3)

Slow Invariant Manifold

The critical manifold is defined by

M 0 := { (⃗ x, ⃗ y) : ⃗ g (⃗ x, ⃗ y, 0) = ⃗ 0 } . ( 4
)
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 a (k -1)-dimensional C ∞ submanifold 5 .
According to Fenichel [START_REF] Fenichel | Persistence and smoothness of invariant manifolds for flows[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] theory if 0 < ε ≪ 1 is sufficiently small, then there exists a function ⃗ G (⃗ y, ε) defined on D such that the manifold

M ε := { (⃗ x, ⃗ y) : ⃗ x = ⃗ G (⃗ y, ε) } , ( 5
)
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 changes6 .

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 [START_REF] Takens | Constrained equations, a study of implicit differential equations and their discontinuous solutions, in Structural stability, the theory of catastrophes and applications in the sciences[END_REF] introduced the "constrained system" defined as follows:

⃗ x = ⃗ f (⃗ x, ⃗ y, 0) , D ⃗ y ⃗ g. ⃗ y = -(D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . (6)
Since, according to Fenichel [START_REF] Fenichel | Persistence and smoothness of invariant manifolds for flows[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], the critical manifold ⃗ g (⃗ x, ⃗ y, 0) may be considered as locally invariant under the flow of system (1), we have:

d⃗ g dt (⃗ x, ⃗ y, 0) = 0 ⇐⇒ D ⃗ x ⃗ g. ⃗ x + D ⃗ y ⃗ g. ⃗ y = ⃗ 0.

By replacing ⃗

x by ⃗ f (⃗ x, ⃗ y, 0) leads to:

D ⃗ x ⃗ g. ⃗ f (⃗ x, ⃗ y, 0) + D ⃗ y ⃗ g. ⃗ y = ⃗ 0.
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:

⃗ x = ⃗ f (⃗ x, ⃗ y, 0) , det(D ⃗ y ⃗ g) ⃗ y = -(adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . (7)

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":

⃗ x = -det(D ⃗ y ⃗ g) ⃗ f (⃗ x, ⃗ y, 0) , ⃗ y = (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) , ⃗ 0 = ⃗ g (⃗ x, ⃗ y, 0) . (8)
where the overdot now denotes the time derivation with respect to τ . Let's notice that José Argémi [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] 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 [START_REF] Benoît | Canards et enlacements[END_REF] as in [START_REF] Benoît | Canards de R 3[END_REF].

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 :

⃗ χ = -det(D ⃗ y ⃗ g) ⃗ f (⃗ χ, ⃗ y, 0) , ⃗ y = (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ χ, ⃗ y, 0) . (9)

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 [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF] for low-dimensional singularly perturbed systems and are defined as singular points of the "normalized slow dynamics" [START_REF] Benoît | Canards et enlacements[END_REF]. 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:

det(D ⃗ y ⃗ g) = 0, (adj(D ⃗ y ⃗ g).D ⃗ x ⃗ g. ⃗ f ) (⃗ x, ⃗ y, 0) = ⃗ 0, ⃗ g (⃗ x, ⃗ y, 0) = ⃗ 0. (10)
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 [START_REF] Argémi | Approche qualitative d'un problème de perturbations singulières dans R 4[END_REF]. 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:

ẋ1 = f 1 (x 1 , x 2 , y 1 , y 2 ) , ( 11a
) ẋ2 = f 2 (x 1 , x 2 , y 1 , y 2 ) , ( 11b
)
ε ẏ1 = g 1 (x 1 , x 2 , y 1 , y 2 ) , ( 11c
)
ε ẏ2 = g 2 (x 1 , x 2 , y 1 , y 2 ) , (11d)
where

⃗ x = (x 1 , x 2 ) t ∈ R 2 , ⃗ y = (y 1 , y 2 ) ∈ R 2 , 0 < ε ≪ 1
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 [START_REF] Chua | Neurons are Poised Near the Edge of Chaos[END_REF] is defined by setting ε = 0 in Eqs. (11c & 11d). Thus, we obtain:

g 1 (x 1 , x 2 , y 1 , y 2 ) = 0, ( 12a
)
g 2 (x 1 , x 2 , y 1 , y 2 ) = 0. ( 12b
)
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:

x 1 = ϕ 1 (x 2 , y 1 , y 2 ) , ( 13a
)
y 1 = ϕ 2 (x 1 , x 2 , y 2 ) . (13b)

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 ):

dg 1 dt = ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 + ∂g 1 ∂y 1 ẏ1 + ∂g 1 ∂y 1 ẏ2 = 0 (14a) dg 2 dt = ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 + ∂g 2 ∂y 1 ẏ1 + ∂g 2 ∂y 1 ẏ2 = 0 (14b)
Eqs. (14a & 14b) may be written as:

∂g 1 ∂y 1 ẏ1 + ∂g 1 ∂y 1 ẏ2 = - ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) (15a) ∂g 2 ∂y 1 ẏ1 + ∂g 2 ∂y 1 ẏ2 = - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2
) (15b)
By solving the system of two equations (15a & 15b) with two unknowns ( ẏ1 , ẏ2 ) we find:

ẏ1 = - ( ∂g1 ∂x1 ẋ1 + ∂g1 ∂x2 ẋ2 ) ∂g2 ∂y2 + ( ∂g2 ∂x1 ẋ1 + ∂g2 ∂x2 ẋ2 ) ∂g1 ∂y2 det [ J (y1,y2) ] , (16a) ẏ2 =
-

( ∂g1 ∂x1 ẋ1 + ∂g1 ∂x2 ẋ2 ) ∂g2 ∂y1 + ( ∂g2 ∂x1 ẋ1 + ∂g2 ∂x2 ẋ2 ) ∂g1 ∂y1 det [ J (y1,y2) ] . ( 16b
)
So, we have the following constrained system:

ẋ1 = f 1 (x 1 , x 2 , y 1 , y 2 ) , ẋ2 = f 2 (x 1 , x 2 , y 1 , y 2 ) , ẏ1 = - ( ∂g1 ∂x1 ẋ1 + ∂g1 ∂x2 ẋ2 ) ∂g2 ∂y2 + ( ∂g2 ∂x1 ẋ1 + ∂g2 ∂x2 ẋ2 ) ∂g1 ∂y2 det [ J (y1,y2) ] , ẏ2 = - ( ∂g1 ∂x1 ẋ1 + ∂g1 ∂x2 ẋ2 ) ∂g2 ∂y1 + ( ∂g2 ∂x1 ẋ1 + ∂g2 ∂x2 ẋ2 ) ∂g1 ∂y1 det [ J (y1,y2) ] , 0 = g 1 (x 1 , x 2 , y 1 , y 2 ) , 0 = g 2 (x 1 , x 2 , y 1 , y 2 ) . ( 17
)

Normalized slow dynamics By rescaling the time by setting

t = -det [ J (y1,y2) ]
τ we obtain the "normalized slow dynamics":

ẋ1 = -f 1 (x 1 , x 2 , y 1 , y 2 ) det [ J (y1,y2) ] = F 1 (x 1 , x 2 , x 3 , y 1 ) , ẋ2 = -f 2 (x 1 , x 2 , y 1 , y 2 ) det [ J (y1,y2) ] = F 2 (x 1 , x 2 , x 3 , y 1 ) , ẏ1 = ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 2 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 2 = G 1 (x 1 , x 2 , x 3 , y 1 ) , ẏ2 = ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 1 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 1 = G 2 (x 1 , x 2 , x 3 , y 1 ) , 0 = g 1 (x 1 , x 2 , y 1 , y 2 ) , 0 = g 2 (x 1 , x 2 , y 1 , y 2 ) . ( 18
)
where the overdot now denotes the time derivation with respect to τ .

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 [START_REF] Guckenheimer | Canards at folded nodes[END_REF] on the tangent bundle of the critical manifold. So, we have:

ẋ2 = -f 2 (x 1 , x 2 , y 1 , y 2 ) det [ J (y1,y2) ] = F 2 (x 2 , y 2 ) , ẏ2 = ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 1 -( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 1 = G 2 (x 2 , y 2 ) . ( 19
)
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:

det [ J (y1,y2) ] = 0, (20a) ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 2 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 2 = 0, (20b) ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 1 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 1 = 0, ( 20c
)
g 1 (x 1 , x 2 , y 1 , y 2 ) = 0, ( 20d
)
g 2 (x 1 , x 2 , y 1 , y 2 ) = 0. ( 20e
)
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" [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] but the "normalized slow dynamics" [START_REF] Guckenheimer | Canards at folded nodes[END_REF].

The Jacobian matrix of system (18) reads:

J (F1,F2,G1,G2) =                ∂F 1 ∂x 1 ∂F 1 ∂x 2 ∂F 1 ∂y 1 ∂F 1 ∂y 2 ∂F 2 ∂x 1 ∂F 2 ∂x 2 ∂F 2 ∂y 1 ∂F 2 ∂y 2 ∂G 1 ∂x 1 ∂G 1 ∂x 2 ∂G 1 ∂y 1 ∂G 1 ∂y 2 ∂G 2 ∂x 1 ∂G 2 ∂x 2 ∂G 2 ∂y 1 ∂G 2 ∂y 2                (21)

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 f 1 (0, 0, 0, 0) = g 1 (0, 0, 0, 0) = 0 ∂g 1 ∂x 2 (0,0,0,0) = ∂g 1 ∂x 3 (0,0,0,0) = ∂g 1 ∂y 1 (0,0,0,0) = 0 [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF] Then, let's make the following assumptions for the non-degeneracy of the folded singularity:

f 2 (0, 0, 0, 0) ̸ = 0 ; ∂g 1 ∂x 1 (0,0,0,0) ̸ = 0 ; ∂ 2 g 1 ∂y 2 1 (0,0,0,0) ̸ = 0. ( 23
)
According to these generic hypotheses Eqs. [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF][START_REF] O'malley | Introduction to Singular Perturbations[END_REF], the Jacobian matrix (21) reads:

J (F1,F2,G1,G2) =          0 0 0 0 -f 2 ∂P ∂x 1 -f 2 ∂P ∂x 2 -f 2 ∂P ∂y 1 -f 2 ∂P
         (24)
where

P = det [ J (y1,y2) ] , a 3i = -f 2 ∂g 2 ∂x 2 ∂ 2 g 1 ∂y 2 ∂x i + ∂g 2 ∂y 2 ( f 2 ∂ 2 g 1 ∂x 2 ∂x i + ∂g 1 ∂x 1 ∂f 1 ∂x i ) for i = 1, 2,
a 3i = -f 2 ∂g 2 ∂x 2 ∂ 2 g 1 ∂y 2 ∂y i + ∂g 2 ∂y 2 ( f 2 ∂ 2 g 1 ∂x 2 ∂y i + ∂g 1 ∂x 1 ∂f 1 ∂y i ) for i = 3, 4,
a 4i = f 2 ∂g 2 ∂x 2 ∂ 2 g 1 ∂y 1 ∂x i - ∂g 2 ∂y 1 ( f 2 ∂ 2 g 1 ∂x 2 ∂x i + ∂g 1 ∂x 1 ∂f 1 ∂x i ) for i = 1, 2,
a 4i = f 2 ∂g 2 ∂x 2 ∂ 2 g 1 ∂y 1 ∂y i - ∂g 2 ∂y 1 ( f 2 ∂ 2 g 1 ∂x 2 ∂y i + ∂g 1 ∂x 1 ∂f 1 ∂y i ) for i = 3, 4.
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:

λ 4 -σ 1 λ 3 + σ 2 λ 2 -σ 3 λ + σ 4 = 0 ( 25
)
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:

λ ( λ 3 -σ 1 λ 2 + σ 2 λ -σ 3 ) = 0 (26)
But, according to Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF], σ 3 vanishes at a pseudo singular point as it's easy to prove it. So, the eigenpolynomial ( 26) is reduced to

λ 2 ( λ 2 -σ 1 λ + σ 2 ) = 0 (27)
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:

σ 1 = T r(J (F1,F2,G1,G2) ) = λ 1 + λ 2 = ∂g 2 ∂x 1 ∂g 1 ∂y 1 ∂f 1 ∂y 2 , σ 2 = 3 ∑ i=1 J ii (F1,F2,G1,G2) = λ 1 λ 2 = ( ∂g 1 ∂y 1 ) 2 [ f 2 2 ( ∂ 2 g 2 ∂x 2 2 ∂ 2 g 2 ∂y 2 2 - ( ∂ 2 g 2 ∂x 2 ∂y 2 ) 2 ) +f 2 ∂g 2 ∂x 1 ( ∂ 2 g 2 ∂y 2 2 ∂f 1 ∂x 2 - ∂ 2 g 2 ∂x 2 ∂y 2 ∂f 1 ∂y 2 )] (28)
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:

C 1 : ∆ = p 2 -4q > 0, C 2 : q < 0. ( 29
)
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 [START_REF] Wechselberger | À propos de canards[END_REF] 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":

ẋ1 = ãx 2 + by 2 + O ( x 1 , ϵ, x 2 2 , x 2 y 2 , y 2 2 ) , ẋ2 = 1 + O (x 1 , x 2 , y 2 , ϵ) , ϵ ẏ1 = cy 1 + O ( ϵx 1 , ϵx 2 , ϵy 2 , x 2 1 , x 2 2 , y 2 2 , x 2 y 2 ) , ϵ ẏ2 = - ( x 1 + y 2 2 ) + O ( ϵx 1 , ϵx 2 , ϵy 2 , ϵ 2 , x 2 1 y 2 , y 3 2 , x 1 x 2 y 2
) .

(

) 30
We establish in Appendix A for any four-dimensional singularly perturbed systems [START_REF] Chua | Neurons are Poised Near the Edge of Chaos[END_REF] with k = 2 slow and m = 2 fast variables that

ã = 1 2 [ f 2 2 ( ∂ 2 g 2 ∂x 2 2 ∂ 2 g 2 ∂y 2 2 - ( ∂ 2 g 2 ∂x 2 ∂y 2 ) 2 ) + f 2 ∂g 2 ∂x 1 ( ∂ 2 g 2 ∂y 2 2 ∂f 1 ∂x 2 - ∂ 2 g 2 ∂x 2 ∂y 2 ∂f 1 ∂y 2 ) ] b = - ∂g 2 ∂x 1 ∂f 1 ∂y 2 , c = ∂g 1 ∂y 1 .
Thus, in his paper Wechselberger [39, p. 3304] provided in the framework of "standard analysis" a generalization of Benoît's theorem [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] 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 [START_REF] Tchizawa | On duck solutions in R 4[END_REF] singular canards perturb to maximal canards solutions for sufficiently small ε ≪ 1.

Proof. See Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF].

Since our method doesn't use the "desingularized vector field" (19) but the "normalized slow dynamics" (18), we have the following proposition:

Proposition 3.
If the normalized slow dynamics [START_REF] Guckenheimer | Canards at folded nodes[END_REF] 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":

ẋ1 = ãx 2 + by 2 + O ( x 1 , ϵ, x 2 2 , x 2 y 2 , y 2 2 ) , ẋ2 = 1 + O (x 1 , x 2 , y 2 , ϵ) , ϵ ẏ1 = cy 1 + O ( ϵx 1 , ϵx 2 , ϵy 2 , x 2 1 , x 2 2 , y 2 2 , x 2 y 2 ) , ϵ ẏ2 = - ( x 1 + y 2 2 ) + O ( ϵx 1 , ϵx 2 , ϵy 2 , ϵ 2 , x 2 1 y 2 , y 3 2 , 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. [START_REF] Tchizawa | On winding duck solutions in R 4[END_REF] it is easy to verify that

σ 1 = T r(J (F1,F2,G1,G2) ) = λ 1 + λ 2 = -bc, σ 2 = 3 ∑ i=1 J ii (F1,F2,G1,G2) = λ 1 λ 2 = 2ãc 2 .
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 [START_REF] Tchizawa | On duck solutions in R 4[END_REF] and Tchizawa [START_REF] Tchizawa | On duck solutions in R 4[END_REF][START_REF] Tchizawa | A direct method for finding ducks in R 4[END_REF][START_REF] Tchizawa | Generic conditions for duck solutions in R 4[END_REF][START_REF] Tchizawa | On the transversality conditions for 4-dim duck solutions[END_REF][START_REF] Tchizawa | On a Local Model for Finding 4-Dimensional Duck Solutions[END_REF][START_REF] Tchizawa | On Relative Stability in 4-Dimensional Duck Solution[END_REF], the existence of "canard solutions" in such system. According to Tchizawa [START_REF] Tchizawa | On the two methods for finding 4-dimensional duck solutions[END_REF], 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 [START_REF] Tchizawa | On the two methods for finding 4-dimensional duck solutions[END_REF] 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. [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] concerning the existence of "canard solutions" in the Hodgkin-Huxley model but with a set of more realistic parameters used in Chua et al. [START_REF] Chua | Hodgkin-Huxley Axon is Made of memristors[END_REF][START_REF] Chua | Neurons are Poised Near the Edge of Chaos[END_REF].

COUPLED FITZHUGH-NAGUMO EQUATIONS

The FitzHugh-Nagumo model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membranes[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF] is a simplified version of the Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] which models in a detailed manner activation and deactivation dynamics of a spiking neuron. By coupling two FitzHugh-Nagumo models Tchizawa & Campbell [START_REF] Tchizawa | On winding duck solutions in R 4[END_REF] and Tchizawa [START_REF] Tchizawa | On duck solutions in R 4[END_REF][START_REF] Tchizawa | On Relative Stability in 4-Dimensional Duck Solution[END_REF] obtained the following four-dimensional singularly perturbed system with two slow and two fast variables:

dx 1 dt = 1 c (y 1 + bx 1 ) , ( 31a
)
dx 2 dt = 1 c (y 2 + bx 2 ) , ( 31b
)
ε dy 1 dt = x 1 - y 3 1 3 + y 2 , ( 31c
)
ε dy 2 dt = x 2 - y 3 2 3 + y 1 . ( 31d
)
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 [START_REF] Tchizawa | A direct method for finding ducks in R 4[END_REF] is defined by setting ε = 0 in Eqs. (31c & 31d). Thus, we obtain:

dx 1 dt = 1 c (y 1 + bx 1 ) , dx 2 dt = 1 c (y 2 + bx 2 ) , dy 1 dt = - 1 c (y 2 + bx 2 ) + y 2 2 c (y 1 + bx 1 ) y 2 1 y 2 2 -1 , dy 2 dt = - 1 c (y 1 + bx 1 ) + y 2 1 c (y 2 + bx 2 ) y 2 1 y 2 2 -1 , 0 = x 1 - y 3 1 3 + y 2 , 0 = x 2 - y 3 2 3 + y 1 .
(32)

Normalized slow dynamics Then, by rescaling the time by setting

t = -det [ J (y1,y2) ] τ = -(y 2 1 y 2 2 -1
) we obtain the "normalized slow dynamics":

dx 1 dt = - 1 c (y 1 + bx 1 ) ( y 2 1 y 2 2 -1 ) = F 1 (x 1 , x 2 , y 1 , y 2 ) , dx 2 dt = - 1 c (y 2 + bx 2 ) ( y 2 1 y 2 2 -1 ) = F 2 (x 1 , x 2 , y 1 , y 2 ) , dy 1 dt = 1 c (y 2 + bx 2 ) + y 2 2 c (y 1 + bx 1 ) = G 1 (x 1 , x 2 , y 1 , y 2 ) , dy 2 dt = 1 c (y 1 + bx 1 ) + y 2 1 c (y 2 + bx 2 ) G 2 (x 1 , x 2 , y 1 , y 2 ) , 0 = x 1 - y 3 1 3 + y 2 , 0 = x 2 - y 3 2 3 + y 1 .
(33)

Pseudo singular points

From Eqs. [START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF], the pseudo-singular points of system (31) are defined by:

det [ J (y1,y2) ] = y 2 1 y 2 2 -1 = 0, (34a) ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 2 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 2 (34b) = 1 c (y 2 + bx 2 ) + y 2 2 c (y 1 + bx 1 ) = 0, ( ∂g 1 ∂x 1 ẋ1 + ∂g 1 ∂x 2 ẋ2 ) ∂g 2 ∂y 1 - ( ∂g 2 ∂x 1 ẋ1 + ∂g 2 ∂x 2 ẋ2 ) ∂g 1 ∂y 1 (34c) = 1 c (y 1 + bx 1 ) + y 2 1 c (y 2 + bx 2 ) = 0, g 1 (x 1 , x 2 , y 1 , y 2 ) = x 1 - y 3 1 3 + y 2 = 0, ( 34d
)
g 2 (x 1 , x 2 , y 1 , y 2 ) = x 2 - y 3 2 3 + y 1 = 0. ( 34e
)
According to Tchizawa & Campbell [START_REF] Tchizawa | On winding duck solutions in R 4[END_REF] and Tchizawa [START_REF] Tchizawa | On duck solutions in R 4[END_REF][START_REF] Tchizawa | A direct method for finding ducks in R 4[END_REF], there are six pseudo singular points, the last four are depending on the parameter b.

(x 1 , x2 , ỹ1 , ỹ2 ) = ( ± 4 3 , ∓ 4 3 , ±1, ∓1 ) , (35a) (x 1 , x2 , ỹ1 , ỹ2 ) =   ± √ 3- √ 9-4b 2 b ( 3 + 2 √ 9 -4b 2 ) 3 √ 2b , ∓ √ 3- √ 9-4b 2 b ( 9 -8b 2 + 3 √ 9 -4b 2 ) 6 √ 2b 2 , ∓ √ 3 - √ 9 -4b 2 2b , ∓ √ 2b √ 3 - √ 9 -4b 2 ) , ( 35b
) (x 1 , x2 , ỹ1 , ỹ2 ) =   ± ( 3 -2 √ 9 -4b 2 ) √ 3+ √ 9-4b 2 b 3 √ 2b , ∓ √ 3+ √ 9-4b 2 b ( 9 -8b 2 -3 √ 9 -4b 2 ) 6 √ 2b 2 , ∓ √ 3 - √ 9 -4b 2 2b , ∓ √ 2b √ 3 - √ 9 -4b 2 ) . (35c)

Canard existence in coupled FitzHugh-Nagumo equations

The Jacobian matrix of system (33) evaluated at the pseudo singular points (35a) reads:

J (F1,F2,G1,G2) =               0 0 2(3 + 4b) 3c - 2(3 + 4b) 3c 0 0 - 2(3 + 4b) 3c 2(3 + 4b) 3c b c b c 1 c -3+8b 3c b c b c - 3 + 8b 3c 1 c               (36)
Remark 4. Although the pseudo singular points have not been shifted at the origin extension of Benoît's generic hypotheses [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF][START_REF] O'malley | Introduction to Singular Perturbations[END_REF] are satisfied. In other words, we have σ 4 = σ 3 = 0.

According to Eqs. ( 28) we find that:

p = σ 1 = T r(J) = + 2 c , q = σ 2 = - 16b(3 + 4b) 9c 2 (37)
Thus, according to Prop. 3, the pseudo singular points are of saddle-type if and only if:

- 16b(3 + 4b) 9c 2 < 0
So, we have the following conditions C 1 and C 2 :

C 1 : ∆ = 4(3 + 8b) 2 9c 2 > 0, C 2 : q = - 16b(3 + 4b) 9c 2 < 0. ( 38
)
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:

b > 0 or b < - 3 4 . ( 39
)
Remark 5. Let's notice that the pseudo singular points are of nodetype if -3 4 < b < 0 as stated by Tchizawa & Campbell [START_REF] Tchizawa | On winding duck solutions in R 4[END_REF] and Tchizawa [START_REF] Tchizawa | On duck solutions in R 4[END_REF][START_REF] Tchizawa | A direct method for finding ducks in R 4[END_REF].

The Jacobian matrix J (F1,F2,G1,G2) of system (33) evaluated at the pseudo singular points (35b) reads:

          0 0 -4 √ 9-4b 2 3c - 2(-9+4b 2 +3 √ 9-4b 2 ) 3bc 0 0 2(-9+4b 2 +3 √ 9-4b 2 ) 3bc 4 √ 9-4b 2 3c 3+ √ 9-4b 2 2c b c 3+ √ 9-4b 2 2bc 3-4 √ 9-4b 2 3c b c 3- √ 9-4b 2 2c 3+4 √ 9-4b 2 3c 3- √ 9-4b 2 2bc           (40)
Remark 6. Although, the pseudo singular points have not been shifted at the origin extension of Benoît's generic hypotheses [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF][START_REF] O'malley | Introduction to Singular Perturbations[END_REF] are satisfied. In other words, we have σ 4 = σ 3 = 0.

According to Eqs. [START_REF] Takens | Constrained equations, a study of implicit differential equations and their discontinuous solutions, in Structural stability, the theory of catastrophes and applications in the sciences[END_REF] we find that:

p = σ 1 = T r(J) = + 3 bc , q = σ 2 = 16 ( 9 -4b 2 ) 9c 2 (41)
Thus, according to Prop. 3, the pseudo singular points are of saddle-type if and only if:

16 ( 9 -4b 2 ) 9c 2 < 0 ∆ = p 2 -4q > 0 and q < 0.
So, we have the following conditions C 1 and C 2 :

C 1 : ∆ = ( 3 bc ) 2 - 64 ( 9 -4b 2 ) 9c 2 > 0, C 2 : q = 16 ( 9 -4b 2 ) 9c 2 < 0. ( 42
)
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:

b > 3 2 or b < - 3 2 . ( 43
)
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 [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] is described by the following system of four nonlinear ordinary differential equations:

dV dt = 1 C M [ I -ḡK n 4 (V -V K ) -ḡNa m 3 h(V -V N a ) -ḡL (V -V L ) ] (44a) dn dt = α n (V )(1 -n) -β n (V )n (44b) dm dt = α m (V )(1 -m) -β m (V )m (44c) dh dt = α h (V )(1 -h) -β h (V )h (44d)
where:

α n (V ) = 0.01(V + 10)/ ( exp V + 10 10 -1 ) , ( 45a
)
β n (V ) = 0.125 exp (V /80) , ( 45b
)
α m (V ) = 0.1(V + 25)/ ( exp V + 25 10 -1 ) , (45c)
β m (V ) = 4 exp(V /18), (45d)
α h (V ) = 0.07 exp(V /20), (45e)
β n (V ) = 1/ ( exp V + 30 10 + 1 ) (45f)
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:

I K = ḡK n 4 (V -V K ) I N a = ḡNa m 3 h(V -V N a ) I L = ḡL (V -V L )
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 , ḡNa 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. (44 & 45) originally chosen by Hodgkin-Huxley and are different from those found in recent literatures where the reference polarity of the voltage V , and the reference direction of the current I are defined as the negative of the voltages and currents. We have opted to adopt the reference assumption in Hodgkin & Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] for ease in comparison of our results with those from Hodgkin and Huxley7 . The parameter values are exactly those chosen in the original Hodgkin-Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] works:

C M = 1.0 µF/cm 2 V N a = -115 mV V K = 12 mV V L = -10.613 mV ḡNa = 120 mS/cm 2 ḡK = 36 mS/cm 2 ḡL = 0.3 mS/cm 2
According to Suckley and Biktashev [START_REF] Suckley | The asymptotic Structure of the Hodgkin-Huxley Equations[END_REF] and Suckley [START_REF] Suckley | The Asymptotic Structure of Excitable Systems[END_REF], dimensionless functions n, h and m 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:

n(V ) = α n (V ) α n (V ) + β n (V ) (46a) h(V ) = α h (V ) α h (V ) + β h (V ) (46b) m(V ) = α m (V ) α m (V ) + β m (V ) (46c) τ n (V ) = 1 α n (V ) + β n (V ) (46d) τ h (V ) = 1 α h (V ) + β h (V ) (46e) τ m (V ) = 1 α m (V ) + β m (V ) (46f)
By using Eqs. 46, the original Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] reads:

dV dt = 1 C M [ I -ḡK n 4 (V -V K ) -ḡNa m 3 h(V -V N a ) -ḡL (V -V L ) ] (47a) dn dt = n -n τ n (47b) dh dt = h -h τ h (47c) dm dt = m -m τ m (47d)
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 [START_REF] Suckley | The asymptotic Structure of the Hodgkin-Huxley Equations[END_REF], Suckley [START_REF] Suckley | The Asymptotic Structure of Excitable Systems[END_REF] and Rubin and Wechselberger [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF], 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 [START_REF] Suckley | The asymptotic Structure of the Hodgkin-Huxley Equations[END_REF], Suckley [START_REF] Suckley | The Asymptotic Structure of Excitable Systems[END_REF] and Rubin and Wechselberger [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] 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 [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] (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. Fig. 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

τ m -1 ≈ 10ms -1 while τ n -1 ≈ τ h -1 ≈ 1ms -1 .
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 [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] 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 ). So, according to Awiszus et al. [START_REF] Awiszus | The singularly perturbed Hodgkin-Huxley equations as a tool for the analysis of repetitive nerve activity[END_REF], Suckley and Biktashev [START_REF] Suckley | The asymptotic Structure of the Hodgkin-Huxley Equations[END_REF], Suckley [START_REF] Suckley | The Asymptotic Structure of Excitable Systems[END_REF] and Rubin and Wechselberger [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] small multiplicative parameters 0 < ε ≪ 1 in the original vector field of the Hodgkin-Huxley Eqs. (47) may be identified while factorizing the right hand side of Eq. (47a) by ḡNa and set: to consistent with the notations of Sec. 3, we obtain:

dx 1 dt = x1 -x 1 τ 1 = f 1 (x 1 , x 2 , y 1 , y 2 ) (48a) dx 2 dt = x2 -x 2 τ 2 = f 1 (x 1 , x 2 , y 1 , y 2 ) (48b) ε dy 1 dt = ȳ1 -y 1 τ 3 = g 1 (x 1 , x 2 , y 1 , y 2 ) (48c) ε dy 2 dt = Ī -ḡK x 4 1 (y 2 -V K ) -ḡNa y 3 1 x 2 (y 2 -V N a ) -ḡL (y 2 -V L ) = g 2 (x 1 , x 2 , y 1 , y 2 ) (48d)
where (x 1 , x2 , ȳ1 ) = (n, h, m) and τ 1,2,3 = τ n,h,m .

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:

x 2 = Ī -ḡK x 4 1 (y 2 -V K ) -ḡL (y 2 -V L ) ḡNa ȳ3 1 (y 2 -V N a ) (49a
)

y 1 = ȳ1 (y 2 ) (49b)
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:

dy 1 dt = ȳ1 -y 1 τ 3 = g 1 (x 1 , x 2 , y 1 , y 2 ) ( 50a
)
dy 2 dt = Ī -ḡK x 4 1 (y 2 -V K ) -ḡNa y 3 1 x 2 (y 2 -V N a ) -ḡL (y 2 -V L ) = g 2 (x 1 , x 2 , y 1 , y 2 ) (50b)
The Jacobian matrix of the fast foliation (50) reads:

J (g1,g2) =    ȳ1 ′ τ 3 -τ ′ 3 ( ȳ1 -y 1 ) τ 2 3 - 1 τ 3 -(ḡ K x 4 + ḡNa y 3 1 x 2 + ḡL ) -3ḡ N a y 2 1 x 2 (y 2 -V N a )    (51)
where the ( ′ ) denotes the derivative with respect to y 2 . Then, taking into account Eqs. (49b), i.e., y 1 = ȳ1 we have:

J (g1,g2) =    ȳ1 ′ τ 3 - 1 τ 3 -(ḡ K x 4 1 + ḡNa ȳ3 1 x 2 + ḡL ) -3ḡ N a ȳ2 1 x 2 (y 2 -V N a )    (52)
So, the determinant of the Jacobian matrix of the fast foliation (50) is:

det ( J (g1,g2) ) = - 1 τ 3 [ ḡK x 4 1 + ḡNa ȳ3 1 x 2 + ḡL + 3ḡ N a ȳ1 ′ ȳ2 1 x 2 (y 2 -V N a ) ] (53)
Thus, the condition for the fold curve is det ( J (g1,g2) ) = 0, which gives:

ḡK x 4 1 + ḡNa ȳ3 1 x 2 + ḡL + 3ḡ N a ȳ1 ′ ȳ2 1 x 2 (y 2 -V N a ) = 0 (54)
Therefore:

x 2 = - ḡK x 4 1 + ḡL ḡNa ȳ2 1 (ȳ 1 + 3 ȳ1 ′ (y 2 -V N a )) (55)
By subtracting Eq. (49a) from Eq. ( 55) we obtain x 1 :

x 1 = x 1f = [ -Ī [ȳ 1 + 3 ȳ1 ′ (y 2 -V N a )] + ḡL (V N a -V L )ȳ 1 + 3 ȳ1 ′ (y 2 -V N a )(y 2 -V L ) ḡK [(V K -V N a )ȳ 1 -3 ȳ1 ′ (y 2 -V N a )(y 2 -V K )] ] 1/4 (56)
Plugging this value of x 1 (56) into Eq. (55) provides:

x 2 = x 2f = Ī + ḡL (V K -V L ) ḡNa ȳ2 1 [(V N a -V K )ȳ 1 + 3 ȳ1 ′ (y 2 -V N a )(y 2 -V K )] (57)
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. [START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF] which reads for the Hodgkin-Huxley model (48): Let's notice that this plot8 is exactly the same as those presented by Rubin and Wechselberger [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF] (Fig. 89) 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 [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF]. Numerical resolution9 of Eq. (60) provides a better approximation of the bifurcation parameter value: ĪC = -4.82988 µA This value corresponds to a voltage y 2 = -3.18136 mV .

For Ī ≈ -4.1, the coordinate of the pseudo singular point can be computed numerically: (x 1 , x 2 , y 1 , y 2 ) = (0.362513, 0.521793, 0.0733782, -2.81908)

According to Proposition 3 we can state that the eigenpolynomial of the Jacobian matrix associated with the "normalized slow dynamics" of the Hodgkin-Huxley model (48) reads:

λ 4 -σ 1 λ 3 + σ 2 λ 2 -σ 3 λ + σ 4 = 0
for which it is easy to prove that σ 4 = σ 3 = 0. So, this eigenpolynomial reduces to:

λ 2 ( λ 2 -σ 1 λ + σ 2 ) = 0
According to Eqs. [START_REF] Takens | Constrained equations, a study of implicit differential equations and their discontinuous solutions, in Structural stability, the theory of catastrophes and applications in the sciences[END_REF] we find that:

p = T r(J) = 144.933, q = σ 2 = -362.924
Thus, according to Prop. 3, the pseudo singular points is of saddletype. Moreover, numerical computation of the eigenvalues of this Jacobian matrix evaluated at the pseudo singular point provides:

(λ 1 , λ 2 , λ 3 , λ 4 ) = (-2.46224, 147.396, 0, 0) So, according to Proposition 3, this pseudo singular point is of saddletype and canard solution may occur in the four-dimensional Hodgkin-Huxley singularly perturbed system (48) for the original set of parameter values.

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 [START_REF] Ginoux | Canards in Memristor's Circuits, Qualitative Theory of Dynamical Systems[END_REF]) 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 [START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF] and then by Szmolyan and Wechselberger [START_REF] Szmolyan | Canards in R 3[END_REF] and finally by Wechselberger [START_REF] Wechselberger | À propos de canards[END_REF] 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.

x = X -Bf 2 2 Y 2 A , y = f 2 Y, z = 1 D [ Z -Cf 2 Y - E G (U -F f 2 Y ) ] , u = U -F f 2 Y G .
(A-4)

The time derivative of system (A-3) gives:

Ẋ = A ẋ + 2By ẏ, Ẏ = ẏ f 2 , Ż = C ẏ + D ż + E u, U = F ż + G u.
(A-5)

Then, multiplying the third and fourth equation of (A-5) by (ε/α) and while replacing in (A-5) ẋ, ẏ, ż and u by the right-hand-side of system (A-2) to:

Ẋ = A ẋ + 2By ẏ, Ẏ = ẏ f 2 , ( ε α ) Ż = ( ε α ) C ẏ + ( ε α ) D ż + ( ε α ) E u, ( ε α ) U = ( ε α ) F ẏ + ( ε α ) G u.
(A-6)

Since ε/α ≪ 1, the first terms of the right-hand-side of the third and fourth equation of (A-16) can be neglected. So we have: ) .

Ẋ = A ( ∂f
(A-7)

Then, by replacing in (A-7) x, y, z and u by the right-hand-side of (A-4) and by identifying with the following system in which we have posed: (ε/α) = ϵ:

Ẋ = ãY + bU + O ( X,

1 FIG. 1 .

 11 FIG. 1.Graph of 1/τ i (ms -1 ) against V (mV ).

  kept as for the original Hodgkin-Huxley model[START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]:C M = 1.0 µF/cm 2 , VNa = -115 mV , VK = 12 mV , VL = -10.613 mV. , h, m, V ) = (x 1 , x 2 , y 1 , y 2 )

2 FIG. 2 .

 22 FIG. 2. Function y 2 ( Ī) for various values of parameter kn = 1, 3, 4.75, 7 exhibiting the the bifurcation parameter value ĪC ≈ -4.8.

FIG. 3 .FIG. 4 .FIG. 5 .

 345 FIG. 3.Phase portrait, canard solution and slow manifold of the Hodgkin-Huxley system (48) in the (n, h, V ) phase space.

  ϵ, Y 2 , Y U, U 2 )This is the result we established in Sec. 2.7.

																					,
			Ẏ = 1 + O (X, Y, U, ϵ) , ϵ Ż = cZ + O ( ϵX, ϵY, ϵU, X 2 , y 2 U, U 2 , Y U ϵ U = -( X + U 2 ) + O ( ϵX, ϵY, ϵU, ϵ 2 , X 2 U, U 3 , XY U ) ,	)	,	(A-8)
	we find:
									ã = A	(	∂f 1 ∂y	-	F G	∂f 1 ∂u	)	f 2 + 2Bf 2 2 ,
									b =	A G	∂f 1 ∂u	,				(A-9)
									c =	∂g 1 ∂z	.
	where
									A = B =	2 1 4	∂g 2 ∂x [ ∂ 2 g 2 ∂ 2 g 2 ∂u 2 , ∂u 2 ∂ 2 g 2 ∂y 2 -	(	∂ 2 g 2 ∂y∂u	) 2 ]	,	(A-10)
									G = -	1 2	∂ 2 g 2 ∂u 2 .
	Finally, we deduce:
	ã =	1 2	[	f 2 2	(	∂ 2 g 2 ∂x 2 2	∂ 2 g 2 ∂y 2 2	-	(	∂ 2 g 2 ∂x 2 ∂y 2	) 2 )	+ f 2	∂g 2 ∂x 1	(	∂ 2 g 2 ∂y 2 2	∂f 1 ∂x 2	-	∂x 2 ∂y 2 ∂ 2 g 2	∂y 2 ∂f 1	) ]
	b = -	∂g 2 ∂x 1	∂f 1 ∂y 2	,
	c =	∂g 1 ∂y 1	.
																					(A-11)

In certain applications these functions will be supposed to be C r , r 1.

It represents the approximation of the slow invariant manifold, with an error of O(ε).

The set D is overflowing invariant with respect to (2) when ε = 0. See Kaper[START_REF] Kaper | An Introduction to Geometric Methods and Dynamical Systems Theory for Singular Perturbation Problems[END_REF] and Jones[START_REF] Jones | Geometric Singular Perturbation Theory in Dynamical Systems[END_REF].

See Benoît[START_REF] Benoît | Systèmes lents-rapides dans R 3 et leurs canards[END_REF][START_REF] Benoît | Perturbation singulière en dimension trois : Canards en un point pseudo singulier noeud[END_REF], Szmolyan and Wechselberger[START_REF] Szmolyan | Canards in R 3[END_REF] and Wechselberger[START_REF] Wechselberger | Existence and Bifurcation of Canards in R 3 in the case of a Folded Node[END_REF][START_REF] Wechselberger | À propos de canards[END_REF].

For more details see Chua et al.[START_REF] Chua | Hodgkin-Huxley Axon is Made of memristors[END_REF][START_REF] Chua | Neurons are Poised Near the Edge of Chaos[END_REF]

The function y 2 ( Ī) solution of (59) has been plotted with Mathematica c ⃝ while using the ContourPlot function used for representing implicit function since such function cannot be expressed explicitly.

This resolution has been made while using the function FindRoot in Mathematica c ⃝ .

ACKNOWLEDGEMENTS

We would like to thank to Ernesto Pérez Chavela for previous discussions related with this work. The authors are partially supported by a MINECO/FEDER grant number MTM2008-03437. The second author is partially supported by a MICINN/FEDER grants numbers MTM2009-03437 and MTM2013-40998-P, by an AGAUR grant number 2014SGR-568, by an ICREA Academia, two FP7+PEOPLE+2012+IRSES numbers 316338 and 318999, and FEDER-UNAB10-4E-378.

allowed this paper to be written. † Partially supported by a DGES grant number PB96-1153.

ȳ1 -y 1 τ 3 = 0, (58a)

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 [START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF], 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.

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.

Normal form of 4D singularly perturbed systems with two fast variables Let's consider the four-dimensional singularly perturbed dynamical system [START_REF] Chua | Neurons are Poised Near the Edge of Chaos[END_REF] with k = 2 slow variables and m = 2 fast and let's make the following change of variables:

By taking into account extension of Benoît's generic hypothesis Eqs. [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF][START_REF] O'malley | Introduction to Singular Perturbations[END_REF] and while using Taylor series expansion the system (11) becomes:

Then, let's make the standard polynomial change of variables: (A-3)

From (A-3) we deduce that: