Bifurcations and canards in the FitzHugh-Nagumo system:
a tutorial in fast-slow dynamics

Bruno F. F. Gonçalves B.F.F. Gonçalves — Centro de Matemática da Universidade do Porto
Rua do Campo Alegre 687, Porto 4169-007, Portugal [email protected] , Isabel S. Labouriau I.S. Labouriau — Centro de Matemática da Universidade do Porto
Rua do Campo Alegre 687, Porto 4169-007, Portugal [email protected] and Alexandre A. P. Rodrigues A.A.P. Rodrigues — Lisbon School of Economics & Management
Centro de Matemática Aplicada à Previsão e Decisão Económica
Rua do Quelhas 6, Lisboa 1200-781, Portugal [email protected]

Abstract.

In this article, we study the FitzHugh-Nagumo $(1,1)$ –fast-slow system where the vector fields associated to the fast/slow equations come from the reduction of the Hodgin-Huxley model for the nerve impulse. After deriving dynamical properties of the singular and regular cases, we perform a bifurcation analysis and we investigate how the parameters (of the affine slow equation) impact the dynamics of the system. The study of codimension one bifurcations and the numerical locus of canards concludes this case-study. All theoretical results are numerically illustrated.

Key words and phrases:

Keywords: FitzHugh-Nagumo system; Fast-slow systems; Canard; Bifurcations.

1. Introduction

Many animals have cells in their various physiological systems, such as the nervous, muscular, and cardiac systems, that are sensitive to certain electrical stimuli. These cells remain at rest most of the time, react to electrical stimulation at a given moment, and then return to a resting state until they are stimulated again. This oscillatory and excitability dynamics is present in neuronal activity and cardiac rhythm – see for instance [19].

In 1952, Hodgkin & Huxley [18] developed a mathematical model that described the propagation of electrical signals along a squid’s giant axon, known as the Hodgkin-Huxley (HH) model. The authors related the excitability of squid nerve cells and the resulting electrical impulse to the potential difference between the inside and outside of the squid’s axon membrane arising from the movement of sodium and potassium ions across the membrane. Using experimental data, they established the propagation of electrical signals along the squid’s giant axon through a nonlinear system of four differential equations [18].

In 1960–1961 Richard FitzHugh [8, 9] created a simplified version of the HH model, considering only two variables. His goal was to simplify the HH model in order to make the dynamics of excitability more perceptible. The system developed by by FitzHugh in 1961 has the following form:

(1)		$\displaystyle\frac{dx}{dt}$	$\displaystyle=a\left(x-\frac{x^{3}}{3}-y+I\right)$
(1)		$\displaystyle\frac{dy}{dt}$	$\displaystyle=\frac{1}{a}\left(x-by-c\right),$

where $a\in{\mathbb{R}}\backslash\{0\},b,c\in{\mathbb{R}}$ are constants and $I$ may be seen as the intensity of a stimulus applied to the axon membrane. Variables $x$ and $y$ of (1) describe the fast evolution of the membrane voltage of a neuron as well as the activation of sodium channels, and the inactivation/activation of the sodium/potassium channels, respectively. They match with the pairs of variables $\left(v,m\right)$ and $\left(h,n\right)$ from the HH model, respectively – cf. [9, 19].

During the sixties, Jin-Ichi Nagumo constructed an electrical circuit that recreated FitzHugh’s system [27]. Nowadays, due to the contributions of both researchers, this model is known as the FitzHugh-Nagumo (FH-N) model. This model has been extensively studied see, for instance, [2] and references therein.

By changing the time scale of system (1), considering $\tau=t/a$ and $\varepsilon=1/a^{2}$ , we obtain an equivalent version of the FH-N model as follows:

(2)		$\displaystyle\varepsilon\frac{dx}{d\tau}$	$\displaystyle=x-\frac{x^{3}}{3}-y+I$
(2)		$\displaystyle\frac{dy}{d\tau}$	$\displaystyle=x-by-c.$

For $0<\varepsilon\ll 1$ , or equivalently, $a\gg 1$ , system (2) is what we call a fast-slow system.

A fast-slow system is a system of differential equations in which some variables have their derivatives with larger magnitude than others. This leads to a system with multiple time scales. The general approach to this type of systems starts by grouping the variables in two disjoint sets: fast variables and slow variables. This separation is introduced in system (2) by the parameter $\varepsilon$ .

In Doss-Bachelet et al. [6] the authors notice that the FH-N model presents two distinct dynamics near a Hopf bifurcation: an attracting focus where all trajectories converge to it if the parameter is less than the bifurcation value; and an attracting relaxation-oscillation periodic solution (limit cycle) if the parameter is greater than the bifurcation point. The work of [6] proceeds with the analysis of the dynamics and applies it to generate bursting oscillations from two coupled FH-N systems, thus leaving the FH-N model little explored.

1.1. Literature

In this section, we review the FH-N equation, as well as its derivations, as one paradigmatic example of mathematical modelling. Although it has been motivated by neuroscience, the FH-N model exists in many variations of the original system and has been studied in several fields due to its simplicity and rich dynamics. Motivated by its broad applicability and generic dynamics, the FH-N system has triggered a huge body of research works characterising their bifurcations and the associated dynamical regimes.

There is an abundance of references in the literature. In what follows, we choose to mention only a few for clarity and the choice is based on our personal preferences. The reader interested in further detail and/or more examples can use the references within those we mention specially those in the review [2].

At the first glance, the simplifications introduced in the passage from the sophisticated HH model to the FH-N model might have a price. Phillipson and Schuster [28] explore numeric and analytically the four-dimensional Hodgkin–Huxley model and demonstrate that, under precise conditions, the FH-N equations can provide quantitative predictions in close agreement with the classic HH equations.

The bifurcations of FH-N system have been extensively analysed as a fast-slow system in [3, 11, 23, 20, 24, 4, 7, 21] uncovering complex patterns of behaviours that emerge from the simple and/or coupled FH-N model. They reveal how changes in the parameters can lead to different dynamical regimes, including equilibria, limit cycles, canards and chaos. Stability techniques allow the identification of critical points where small perturbations may lead to qualitatively different behaviours, shedding light on the robustness and sensitivity of excitable systems. Systems exhibiting excitability often display a threshold behaviour, where a certain level of input of perturbation is required to trigger the response – cf. [6].

We would like to refer the paper [10] for the study, not using the fast-slow approach, of (1), where global bifurcations are presented in the context of the global bifurcation diagram. Codimension one and two bifurcations of Hopf, homoclinic, saddle-node and Bogdanov-Takens are obtained.

The work [1] investigates the complex phenomena of canard explosion with mixed-mode oscillations, which can be observed in a fractional-order FH-N model. The study has highlighted the appearance of patterns of solutions with an increasing number of small-amplitude oscillations in each of such patterns, as one parameter of the fractional-order FH-N model is varied.

The paper [25] extends the FH-N model by including a time-varying threshold between electrical “silence” and electrical firing, as well as a periodic forcing for the intensity of the stimulus. Authors provide sufficient conditions for the existence of periodic solutions in such differential systems. This may be related with the work [31], where the author explores the dynamical behaviour of non-autonomous, almost-periodic discrete FH-N system defined on infinite lattices. The non-autonomous infinite-dimensional system has a uniform attractor which attracts all solutions uniformly with respect to the translations of external terms. Authors establish the upper semi-continuity of uniform attractors when the infinite-dimensional system is approached by a family of finite-dimensional systems.

In [5] authors propose a new method for the identification of the FH-N model dynamics via deterministic learning (with diffusion). The FH-N model is then described by partial differential equations which are transformed into a set of ordinary differential equations whose dynamics has been studied. Authors found spiral waves and an accurate identification of dynamics underlying the recurrent trajectory corresponding to any spatial point is achieved.

Spiking and bursting patterns observed in nerve membranes seem to be important when we investigate information representation model in the brain. Many topologically different bursting responses are observed in the mathematical models and their related bifurcation mechanisms have been clarified in [30]. In this article, the authors propose a design method to generate bursting responses in FitzHugh–Nagumo model with a simple periodic external force based on bifurcation analysis. Some effective parameter perturbations for the amplitude of the external input are given from the two-parameter bifurcation diagram.

1.2. Novelty and structure of the article

Our goal with the present article is to exhibit the various dynamics in the multiple time scale model and motivate their existence in the light of methods from geometric singular perturbation theory and bifurcation theory. This article studies the dynamics of a fast-slow system derived from the FH-N model according to the geometric singular perturbation theory and the bifurcation theory methods. We provide an analytic proof that the fast-slow FH-N system presents relaxation oscillation dynamics as well as periodic solutions induced by Hopf bifurcation – we fully analyse the type of bifurcations and we check rigorously all conditions for the emergence of a homoclinic orbit and canards connecting these two distinct dynamics. We illustrate each result with numerical computations.

This article is organised as follows. In Section 3 we introduce general concepts of fast-slow systems and geometric singular perturbation theory and apply them to the analysis of the singular case of the FH-N system. A relaxation-oscillation periodic solution is studied with an estimate of its period in Section 3.2. In Section 4 we describe Fenichel’s theorem and link singular to regular dynamics and the way the latter converges to the former as $\varepsilon\to 0$ applying it to FH-N. In Section 5 we perform a bifurcation analysis of FH-N and investigate the influence of parameters in the qualitatively different dynamics of the model. The study of singularities and the emergence of canard trajectories in Section 6 concludes this tutorial. Section 7 discusses the results presented and states the natural continuation of the present work.

We have endeavoured to make a self contained exposition bringing together all topics related to the proofs. We have drawn illustrative figures to make the paper easily readable. All figures in this article were created through numerical simulations conducted in Matlab, using integration functions such as ode15s or ode23s, and, in more sensitive examples, the MatCont toolbox developed by Willy Govaerts, Yuri A. Kuznetsov and Hil G.E. Meijer [14].

2. Multiple Time Scale Theory and the FitzHugh-Nagumo System

In this section, general concepts of two time scale systems are introduced and illustrated with the FH-N system. Further results on fast-slow systems may be found in [24]. We start by introducing some necessary definitions.

Definition 2.1.

For $m,n\in{\mathbb{N}}$ , a $(m,n)$ –fast-slow system is a system of $m+n$ ordinary differential equations of the form:

(3)			$\displaystyle\varepsilon\,\dfrac{\mathrm{d}x}{\mathrm{d}\tau}=\varepsilon\dot{x}=f(x,y,\varepsilon),$
(3)			$\displaystyle\phantom{\varepsilon\,}\dfrac{\mathrm{d}y}{\mathrm{d}\tau}=\phantom{\varepsilon}\dot{y}=g(x,y,\varepsilon),$

where $f\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{m}\times{\mathbb{R}}^{n}\times{\mathbb{R}}\to{\mathbb{R}}^{m}$ , $g\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{m}\times{\mathbb{R}}^{n}\times{\mathbb{R}}\to{\mathbb{R}}^{n}$ are $C^{2}$ and $0<\varepsilon\ll 1.$ We denote the fast variable(s) by $x\in{\mathbb{R}}^{m}$ and $y\in{\mathbb{R}}^{n}$ is called the slow variable(s). Considering $t=\tau/\varepsilon$ , system (3) is equivalent to:

(4)			$\displaystyle\dfrac{\mathrm{d}x}{\mathrm{d}t}=x^{\prime}=\phantom{\varepsilon\,}f(x,y,\varepsilon),$
(4)			$\displaystyle\dfrac{\mathrm{d}y}{\mathrm{d}t}=y^{\prime}=\varepsilon\,g(x,y,\varepsilon).$

The fast time scale is represented by $t$ while $\tau$ refers to the slow time scale. In order to better understand the concepts and the terminology associated with these systems, the following definitions and results are illustrated by the FH-N (1,1)–fast-slow system.

Example 2.1.

As described in [6], consider the FH-N (1,1)–fast-slow system:

(5)	$\displaystyle\text{Slow time scale }\left(\tau\right)\qquad\qquad$	$\displaystyle\text{Fast time scale }{\textstyle\left(t=\tau/\varepsilon\right)}$
	$\displaystyle\varepsilon\,\dfrac{\mathrm{d}x}{\mathrm{d}\tau}=\varepsilon\dot{x}=f(x,y,\varepsilon)\qquad\qquad$	$\displaystyle\dfrac{\mathrm{d}x}{\mathrm{d}t}=x^{\prime}=\phantom{\varepsilon}\,f(x,y,\varepsilon)$
	$\displaystyle\phantom{\varepsilon\,}\dfrac{\mathrm{d}y}{\mathrm{d}\tau}=\phantom{\varepsilon}\dot{y}=g(x,y,\varepsilon)\qquad\qquad$	$\displaystyle\dfrac{\mathrm{d}y}{\mathrm{d}t}=y^{\prime}=\varepsilon\,g(x,y,\varepsilon),$

where $f,g\mathrel{\mathop{\mathchar 58\relax}}\;{\mathbb{R}}^{3}\longrightarrow{\mathbb{R}}$ are given by

f(x,y,\varepsilon)=-y+4x-x^{3},\quad g(x,y,\varepsilon)=x-by-c

with parameters $b,c\in{{\mathbb{R}}}$ and $0<\varepsilon\ll 1$ . $\diamond$

Both cubic and affine functions of Example 2.1, $f$ and $g$ , come naturally from the reduction of the HH model [19]. Note that in FH-N the functions $f$ and $g$ do not depend on $\varepsilon$ . System (LABEL:fhn) will be referred to as the FH-N system in two variables.

3. The Singular Case

When analysing a fast-slow system, it is useful to consider the singular case $\varepsilon=0$ . Techniques such as nullclines analysis provide a geometric understanding of the model’s behaviour.

Definition 3.1.

For $\varepsilon=0$ , system (3) in the slow time scale becomes a system of differential-algebraic equations:

(6)			$\displaystyle 0=f(x,y,\varepsilon),$
(6)			$\displaystyle\dot{y}=g(x,y,\varepsilon).$

System (6) is called slow subsystem, its equations are referred to as the reduced equations and its flow the slow/reduced flow.

Definition 3.2.

For $\varepsilon=0$ , system (4) in fast time is called fast subsystem:

(7)			$\displaystyle x^{\prime}=f(x,y,\varepsilon),$
(7)			$\displaystyle y^{\prime}=0.$

The set of equations in system (7) is referred to as the layer equations and its flow as the fast flow.

Systems (6) and (7) allow the understanding of the dynamics of (3) and (4) for small $\varepsilon$ , $0<\varepsilon\ll 1$ , since the flow may be seen as either the fast or the slow flow, depending on the region in the phase space. Near the set defined by $f\left(x,y,0\right)=0$ , one expects the flow to be $C^{2}$ –close to the slow flow as we will see later in Section 4.

Definition 3.3.

For $m,n\in{\mathbb{N}}$ , given a $(m,n)$ –fast-slow system as in (2.1), the critical manifold is the set

(8)

{}\mathcal{C}_{0}=\left\{(x,y)\in{\mathbb{R}}^{m}\times{\mathbb{R}}^{n}\mathrel{\mathop{\mathchar 58\relax}}f(x,y,0)=0\right\}.

Example 3.1.

The critical manifold for the FH-N system (LABEL:fhn) is

\mathcal{C}_{0}=\left\{(x,y)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}f(x,y,0)=0\right\}=\left\{(x,y)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}y=4x-x^{3}\right\},

a cubic shaped curve, parametrised by

(9)

x\mapsto(x,4x-x^{3})=(x,y),

as depicted in Figure (1(a)). The map $\varphi(x)=4x-x^{3}$ is odd $\left(\forall x\in{\mathbb{R}},-\varphi\left(x\right)=\varphi\left(-x\right)\right)$ . $\diamond$

All points in $\mathcal{C}_{0}$ are equilibria of the fast subsystem. Through Hartman-Grobman’s Theorem [17], we may define the (normal) hyperbolicity and stability of the points in $\mathcal{C}_{0}$ with respect to the fast variables, in the following way:

Definition 3.4.

For $m,n\in{\mathbb{N}}$ , given a $(m,n)$ –fast-slow system with critical manifold $\mathcal{C}_{0}$ :

(1)

The subset $S\subset\mathcal{C}_{0}$ is said to be normally hyperbolic if for all $\mathfrak{p}\in{S}$ , we have

$\left(\mathsf{D}_{x}f\right)\left(\mathfrak{p},0\right)\text{ has no eigenvalues with zero real part,}$

where $\left(\mathsf{D}_{x}f\right)$ denotes the Jacobian matrix of $f$ with respect to $x$ .
(2)

A normally hyperbolic set $S$ is said to be attracting/repelling if, for all $\mathfrak{p}\in{S}$ , all eigenvalues of $\left(\mathsf{D}_{x}f\right)\left(\mathfrak{p},0\right)$ have negative/positive real part, respectively. It is of saddle type if $\left(\mathsf{D}_{x}f\right)\left(\mathfrak{p},0\right)$ has eigenvalues with both positive and negative real part.

When considering a (1,1)–fast-slow system, since there is only one fast variable, we have $\mathsf{D}_{x}f=\dfrac{\partial f}{\partial x}$ . For convenience, we use the notation $\varphi_{\sigma}$ when referring to the partial derivative of a real map $\varphi$ in order to its variable $\sigma$ .

Definition 3.5.

Given a (1,1)–fast-slow system with critical manifold $\mathcal{C}_{0}$ , the point $\mathfrak{p}\in\mathcal{C}_{0}$ is said to be a fold point if:

f_{x}\left(\mathfrak{p},0\right)=0\,,\quad f_{xx}\left(\mathfrak{p},0\right)\neq 0\quad\text{ and }\quad f_{y}\left(\mathfrak{p},0\right)\neq 0.

If $g\left(\mathfrak{p},0\right)\neq 0$ , the fold point is called regular.

Example 3.2.

With respect to (LABEL:fhn), we have $f\left(x,y,0\right)=-y+4x-x^{3}$ and, for $\mathfrak{p}\in{\mathcal{C}_{0}}$ , one has:

f_{x}\left(\mathfrak{p},0\right)=4-3x^{2}\neq 0\qquad\text{ if }\,\,x\neq\pm\dfrac{2}{\sqrt{3}}.

Therefore, the equilibria $P_{+}=\left(2/\sqrt{3},\,16/3\sqrt{3}\right)$ and the set $P_{-}=\left(-2/\sqrt{3},\,-16/3\sqrt{3}\right)$ are fold points and $S_{0}=\mathcal{C}_{0}\setminus\left\{\left(x,\varphi\left(x\right)\right)\mathrel{\mathop{\mathchar 58\relax}}\,x=\pm 2/\sqrt{3}\right\}$ is normally hyperbolic. We may split $S_{0}$ into three distinct open subsets:

	$\displaystyle\mathcal{C}_{0\text{L}}$	$\displaystyle=\mathcal{C}_{0}\cap{\textstyle\left\{\left(x,y\right)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}x<\dfrac{-2}{\sqrt{3}}\right\}},$
	$\displaystyle\mathcal{C}_{0\text{M}}$	$\displaystyle=\mathcal{C}_{0}\cap{\textstyle\left\{\left(x,y\right)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}\dfrac{-2}{\sqrt{3}}<x<\dfrac{2}{\sqrt{3}}\right\}}\quad\text{and}$
	$\displaystyle\mathcal{C}_{0\text{R}}$	$\displaystyle=\mathcal{C}_{0}\cap{\textstyle\left\{\left(x,y\right)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}x>\dfrac{2}{\sqrt{3}}\right\}}.$

Therefore, we get

(10)

S_{0}=\mathcal{C}_{0\text{L}}\cup\mathcal{C}_{0\text{M}}\cup\mathcal{C}_{0\text{R}},

where $\mathcal{C}_{0\text{L}}$ and $\mathcal{C}_{0\text{R}}$ are attracting subsets and $\mathcal{C}_{0\text{M}}$ is repelling. This means that in the fast flow, trajectories of (LABEL:fhn) move towards either $\mathcal{C}_{0\text{L}}$ or $\mathcal{C}_{0\text{R}}$ . $\diamond$

Since the critical manifold $\mathcal{C}_{0}$ is parametrised by the fast variable $x$ (cf. Eq.(9)), it is possible to explicitly define the slow flow through the fast variable $x$ . This allows the stability analysis of the fast and slow flows expressed in terms of the fast variable.

In system (LABEL:fhn), by differentiating $f(x,y,0)=-y+4x-x^{3}=0$ with respect to $\tau$ we have:

\dfrac{\partial f}{\partial x}\dot{x}+\dfrac{\partial f}{\partial y}\dot{y}=(4-3x^{2})\dot{x}-\dot{y}=0.

Since $\dot{y}=x-by-c$ , then

\dot{x}=\dfrac{x-by-c}{(4-3x^{2})}.

Taking into account that $y=4x-x^{3}$ , we may conclude that

(11)

\dot{x}=\dfrac{bx^{3}+(1-4b)x+c}{4-3x^{2}}\eqqcolon\psi\left(x,b,c\right).

As expected, this explicit form of the slow subsystem is not defined at the fold points $x=\pm 2/\sqrt{3}$ .

If the fold points are regular (cf. Def. 3.5), then all equilibria $\left(x^{*},y^{*}\right)\in\mathcal{C}_{0}$ of the system (LABEL:fhn) belong necessarily to $S_{0}$ . Their stability in the slow direction may be determined in the following way:

\dfrac{\partial\psi}{\partial x}=\psi_{x}\left(x,b,c\right)=\dfrac{1-b\left(4-3x^{2}\right)+6x\,\psi\left(x,b,c\right)}{\left(4-3x^{2}\right)}.

Since $\psi\left(x^{*},b,c\right)=0$ , then

\psi_{x}\left(x^{*},b,c\right)=\dfrac{1-b\left(4-3{x^{*}}^{2}\right)}{4-3{x^{*}}^{2}}\,.

For $b>0$ , if $|x^{*}|>2/\sqrt{3}$ , the equilibrium attracts the slow flow. Otherwise, if $|x^{*}|<2/\sqrt{3}$ the equilibrium may either attract or repel the slow flow depending on the value of $1-b\left(4-3{x^{*}}^{2}\right)$ being negative or positive, respectively. Thus, the equilibrium point $\left(x^{*},y^{*}\right)$ of (LABEL:fhn) is stable if $\left(x^{*},y^{*}\right)\in\mathcal{C}_{0\text{L}}\cup\mathcal{C}_{0\text{R}}$ and unstable if $\left(x^{*},y^{*}\right)\in\mathcal{C}_{0\text{M}}$ .

For $b<0$ , the stability in the slow direction is the opposite of the previous one. Hence, the equilibrium point $\left(x^{*},y^{*}\right)$ is unstable if $\left(x^{*},y^{*}\right)\in\mathcal{C}_{0\text{M}}$ and either stable or unstable of saddle type if $\left(x^{*},y^{*}\right)\in\mathcal{C}_{0\text{L}}\cup\mathcal{C}_{0\text{R}}$ .

In Example 3.2, we have disregarded the case where a fold point is also an equilibrium of the system, i.e., $x^{*}=\pm 2/\sqrt{3}$ . This particular situation requires special attention and will be analysed later in Section 6. For the remainder of this section, we describe the dynamics in the cases where all equilibria lie in the set $S_{0}$ defined in (10).

3.1. Dynamics

Depending on the values of $b$ and $c$ , there exists at least one and at most three equilibria, as result of the intersection of the cubic curve, $f\left(x,y,0\right)=0$ , with the line, $\dot{y}=0$ (cf. Figure 1(b)). Moreover, when the system has three equilibria only the following scenarios may occur: either all three equilibria belong to $\mathcal{C}_{0\text{M}}$ or else they are each one in a distinct region, $\mathcal{C}_{0\text{L}}$ , $\mathcal{C}_{0\text{M}}$ and $\mathcal{C}_{0\text{R}}$ .

Fix $A=\left(x_{a},y_{a}\right)\in{\mathbb{R}}^{2}\setminus\mathcal{C}_{0}$ . Consider the half trajectory $\gamma_{A}$ for $t\geq 0$ that starts at $A$ . Since $y^{\prime}=0$ , the fast flow runs horizontally towards $\mathcal{C}_{0\text{L}}$ or $\mathcal{C}_{0\text{R}}$ . Without loss of generality, suppose that $\gamma_{A}$ moves towards $\mathcal{C}_{0\text{L}}$ and let $B=\left(x_{b},y_{a}\right)\in\mathcal{C}_{0\text{L}}$ be the intersection point of $\gamma_{A}$ with $\mathcal{C}_{0}$ . Here, five different scenarios may occur (cf. Figure 2):

(i)

The point $B$ is an equilibrium;
(ii)

The point $B$ is not an equilibrium, but there exists a stable equilibrium $C=\left(x_{c},y_{c}\right)\in\mathcal{C}_{0\text{L}}\setminus{B}$ ;
(iii)

The point $B$ is not an equilibrium, but there exists an unstable equilibrium $C=\left(x_{c},y_{c}\right)\in\mathcal{C}_{0\text{L}}\setminus{B}$ ;
(iv)

There is no equilibrium in $\mathcal{C}_{0\text{L}}$ , but it exists in $\mathcal{C}_{0\text{R}}$ ;
(v)

All equilibria lie in $\mathcal{C}_{0\text{M}}$ .

Scenario (i) is trivial. If $B=\left(x_{b},y_{a}\right)$ is an equilibrium of (LABEL:fhn), then by definition $\left(\dot{x},\dot{y}\right)=0$ , and therefore, the trajectory initiated at $A$ accumulates at $B$ (cf. Figure 2(a)).

Scenario (ii) is not much different. Upon reaching the point $B$ , the trajectory $\gamma_{A}$ moves through the slow flow of the system defined by $\dot{y}=x-by-c$ . Since the equilibrium point $C$ is attracting, and since both points $B$ and $C$ are in $\mathcal{C}_{0\text{L}}$ , the trajectory $\gamma_{A}$ moves from $B$ to $C$ along the branch $\mathcal{C}_{0\text{L}}$ of the critical manifold, where it accumulates (cf. Figure 2(b)).

Scenario (iii): First note that $\dfrac{\partial g}{\partial y}\left(x,y,0\right)=-b>0$ is constant for all $(x,y)\in\mathcal{C}_{0}$ . Since $\mathcal{C}_{0\text{L}}$ and $\mathcal{C}_{0\text{R}}$ attract the slow flow we have that any equilibrium point in either $\mathcal{C}_{0\text{L}}$ or $\mathcal{C}_{0\text{R}}$ is a saddle. In particular, the equilibrium point $C$ is a saddle.

Since $\mathcal{C}_{0\text{L}}$ attracts the fast flow, we know that $\gamma_{A}$ follows $\mathcal{C}_{0\text{L}}$ away from $C$ . Thus, if $y_{a}>y_{c}$ , then $\gamma_{A}$ remains in $\mathcal{C}_{0\text{L}}$ and the second coordinate of $\gamma_{A}\left(\tau\right)$ tends to $+\infty$ as in Figure 2(c).

However, if $y_{a}<y_{c}$ , the trajectory follows the slow flow to the fold point $P_{-}=\left(-2/\sqrt{3},\,-16/3\sqrt{3}\right)$ . Here, since the slow flow is not defined at the fold points (cf. (11)), the trajectory $\gamma_{A}$ moves through the fast flow, horizontally, to the point $D=\left(4/\sqrt{3},\,-16/3\sqrt{3}\right)\in\mathcal{C}_{0\text{R}}$ as in Figure 2(d). In this case there are two possibilities.

The first possibility is that there is an equilibrium point $E=\left(x_{e},y_{e}\right)\in\mathcal{C}_{0\text{R}}$ and that both $y_{c}$ and $y_{e}$ lie outside the interval $\left(-16/3\sqrt{3},\,\,16/3\sqrt{3}\right)$ delimited by the second coordinates of the two fold points. Then $\gamma_{A}$ follows the slow flow up to the other fold point $P_{+}=\left(2/\sqrt{3},\,\,16/3\sqrt{3}\right)$ , where it jumps back to $\mathcal{C}_{0\text{L}}$ and continues cycling around between the two branches. The trajectory $\gamma_{A}$ presents dynamics similar to case (v).

In any other case, i.e., if either there is no equilibrium in $\mathcal{C}_{0\text{R}}$ or if one of the equilibria has its second coordinate in the interval $\left(-16/3\sqrt{3},\,\,16/3\sqrt{3}\right)$ , then the second coordinate of $\gamma_{A}\left(\tau\right)$ tends to $\pm\infty$ as $\tau$ tends to $+\infty$ . This is because, if $E=\left(x_{e},y_{e}\right)\in\mathcal{C}_{0\text{R}}$ with $y_{e}<-16/3\sqrt{3}$ then the trajectory follows $\mathcal{C}_{0\text{R}}$ up to the fold point $P_{+}$ , jumps to $\mathcal{C}_{0\text{L}}$ and follows it up, away from $C$ and with its second coordinate going to $+\infty$ . If either there is no equilibrium in $\mathcal{C}_{0\text{R}}$ or if $E=\left(x_{e},y_{e}\right)\in\mathcal{C}_{0\text{R}}$ with $y_{e}>-16/3\sqrt{3}$ , then $\gamma_{A}$ follows $\mathcal{C}_{0\text{R}}$ down with its second coordinate going to $-\infty$ as in Figure 2(d).

Scenario (iv): let $E=\left(x_{e},y_{e}\right)$ be an equilibrium point in $\mathcal{C}_{0\text{R}}$ . Geometrically, the point $E$ is the intersection of the cubic curve $x^{\prime}=0$ with the line $\dot{y}=0$ . Since there are no more equilibria in $\mathcal{C}_{0\text{L}}$ , we may observe that, if $E$ is unstable, then $y_{a}>y_{e}$ . The proof is simple: if $E$ is unstable, then $\psi_{x}<0$ where $\psi$ was defined in (11). Since $f_{x}\left(E,0\right)>0$ , it follows that $b<0$ and therefore the line has a negative slope. Now, assume by contradiction that $y_{b}<y_{e}$ . The line $\dot{y}=0$ necessarily intersects the cubic surface in $\mathcal{C}_{0\text{L}}$ and hence, there exists an equilibrium in $\mathcal{C}_{0\text{L}}$ giving our contradiction. Thus, at $B$ , we have $\dot{y}>0$ and therefore, the trajectory $\gamma_{A}$ remains in $\mathcal{C}_{0\text{L}}$ , with the first coordinate of $\gamma_{A}\left(\tau\right)$ going to $-\infty$ as $\tau\rightarrow\infty$ while the second coordinate of $\gamma_{A}\left(\tau\right)$ goes to $\infty$ as in Figure 2(f). If the point $E$ is stable, $\gamma_{A}$ travels along $\mathcal{C}_{0\text{L}}$ from $B$ to $P_{-}$ , jumps in fast time to point $D$ in the branch $\mathcal{C}_{0\text{R}}$ , ending at $E$ (cf. Figure 2(e)).

Scenario (v): by using the symmetry of $f$ ( $f$ is odd), one may check easily that the parameter $b$ is positive and that there are either three equilibria or only one in $\mathcal{C}_{0\text{M}}$ . In both cases, since the points belong to the unstable region of $\mathcal{C}_{0}$ , the dynamics of $\gamma_{A}$ does depend on the number of equilibria in the system. Therefore, we can observe that any point in the region $\mathcal{C}_{0\text{L}}$ satisfies $\dot{y}<0$ , so the trajectory $\gamma_{A}$ moves along $\mathcal{C}_{0\text{L}}$ from point $B$ to the fold point $P_{-}$ , and then jumps, in fast time, to point $D$ . Similarly, $\gamma_{A}$ continues to $P_{+}$ and then jumps to point $F$ . Thus, $\gamma_{A}$ goes through points $F,P_{-},D$ and $P_{+}$ periodically and, therefore, $\gamma_{A}$ is a periodic solution. The trajectory neither diverges nor converges to an equilibrium point (cf. Figure 3(a)).

Any trajectory starting at any point in ${\mathbb{R}}^{2}\setminus\mathcal{C}_{0}$ moves in fast time to a point in one of the stable regions of the critical manifold $\mathcal{C}_{0}$ . In the previous scenarios, we have assumed that this point belongs to the branch $\mathcal{C}_{0\text{L}}$ , but the dynamics for the case where the point is located on the branch $\mathcal{C}_{0\text{R}}$ is analogous to the former.

3.2. Period of the limit cycle

In Scenario (iii) when the equilibria of $\mathcal{C}_{0\text{L}}\cup\mathcal{C}_{0\text{R}}$ satisfy $|y|>\frac{16}{3\sqrt{3}}$ and in Scenario (v), we have observed the emergence of a periodic solution alternating between the slow and the fast flows. Now, in order to estimate its period, we are going to determine the points responsible for linking slow and fast flows.

The transition from the slow flow to the fast one occurs at the fold points $P_{\pm}=\pm\textstyle\left(2/\sqrt{3},\,16/3\sqrt{3}\right)$ . Similarly, the transition from the fast flow to the slow one happens at the points ${F=\textstyle\left(-4/\sqrt{3},\,16/3\sqrt{3}\right)}$ and ${\textstyle D=\left(4/\sqrt{3},\,-16/3\sqrt{3}\right)}$ .

Let $\gamma_{0}$ be the periodic trajectory of system (LABEL:fhn) with initial condition at $D$ (cf. Figure 3). As we have seen, $\gamma_{0}$ follows periodically the cycle:

(12)

\gamma_{0}\mathrel{\mathop{\mathchar 58\relax}}D\;\overset{\text{slow flow}}{\longrightarrow}\;P_{+}\;\overset{\text{fast flow}}{\longrightarrow}\;F\;\overset{\text{slow flow}}{\longrightarrow}\;P_{-}\;\overset{\text{fast flow}}{\longrightarrow}\;D\;\overset{\text{slow flow}}{\longrightarrow}\;\left(...\right)\;.

We intend to compute the period $\mathcal{T}_{\gamma_{0}}$ of this cycle. Since $\varepsilon=0$ , then the fast flow of the system is traversed almost instantaneously. Therefore, $\mathcal{T}_{\gamma_{0}}$ is determined by the durations $\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{L}}}$ and $\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}$ of $\gamma_{0}$ along the branches $\mathcal{C}_{0\text{L}}$ and $\mathcal{C}_{0\text{R}}$ , respectively. Since $f$ is odd, we have $\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{L}}}=\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}$ . Thus,

(13)

{}\mathcal{T}_{\gamma_{0}}\approx\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{L}}}+\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}=2\;\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}.

Hence, the time spent by the trajectory in $\mathcal{C}_{0\text{R}}$ is given by:

\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}\equiv\int_{0}^{\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}}\;\mathrm{d}\tau.

Using (11) we determine $\mathcal{T}_{\mathcal{C}_{0\text{R}}}$ with respect to the fast variable:

(14)

\mathcal{T}_{\gamma_{0}}\approx 2\;\mathcal{T}_{\gamma_{0}}^{\mathcal{C}_{0\text{R}}}=2\int_{\frac{4}{\sqrt{3}}}^{\frac{2}{\sqrt{3}}}\dfrac{4-3x^{2}}{bx^{3}+(1-4b)x+c}\;\mathrm{d}x,

that may be evaluated for any suitable value of $b$ and $c$ . For instance, for $b=c=0$ we are in case (v) and

\mathcal{T}_{\gamma_{0}}=2\int_{\frac{4}{\sqrt{3}}}^{\frac{2}{\sqrt{3}}}\frac{4-3x^{2}}{x}\;\mathrm{d}x=12-8\log\left(2\right)\approx 6.45\,.

In this section we have described the dynamics of (LABEL:fhn) in the singular case $\varepsilon=0$ , in the various cases according to the number and stability of its equilibria, as shown in Figures 2 and 3. We have also obtained conditions under which (LABEL:fhn) with $\varepsilon=0$ has a periodic solution that alternates between the fast and the slow flows and we have determined its period. In the next section we look at the dynamics of (3) in the nonsingular case, i.e., when $\varepsilon\overset{>}{\sim}0$ .

4. The Regular Case

The independence of the time scales when $\varepsilon=0$ makes the system simpler to understand. By allowing the dynamics to be separated into two distinct phases, slow and fast, we can study each one independently and understand its underlying dynamics.

As we will see in Fenichel’s theorem (Theorem 1 stated below) when $0<\varepsilon\ll 1$ , the dynamics of the perturbed system is similar to that of the singular system, with a deviation of $\mathcal{O}\left(\varepsilon\right)$ , where $\mathcal{O}$ stands for the usual Landau notation. That is, the smaller the $\varepsilon$ , the more similar the system trajectories are to those described in the singular system (cf. Figure 4(a)). Fenichel’s theorem guarantees that for $0<\varepsilon\ll 1$ , there exists a set $S_{\varepsilon}$ very close to $S_{0}$ that exhibits the same behaviour as $S_{0}$ .

Theorem 1 (Fenichel, 1979 [7]).

Suppose $S_{0}$ is a compact normally hyperbolic submanifold of the critical manifold $\mathcal{C}_{0}$ of (3) and that $f,g\in{C^{r}}\,\,\left(2\leq r<\infty\right)$ . Then for $\varepsilon>0$ , sufficiently small, the following hold:

(H1)

There exists a locally invariant manifold $S_{\varepsilon}$ diffeomorphic to $S_{0}$ .
(H2)

The set $S_{\varepsilon}$ has Hausdorff¹¹1The Hausdorff distance between two nonempty sets $V,W\subset{\mathbb{R}}^{m+n}$ is defined by $d_{H}\left(V,W\right)\mathrel{\mathop{\mathchar 58\relax}}=\max\left\{\sup_{v\in V}\inf_{w\in W}||v-w||,\,\sup_{w\in W}\inf_{v\in V}||v-w||\right\}.$ distance $\mathcal{O}\left(\varepsilon\right)$ , as $\varepsilon\to 0$ , from $S_{0}$ .
(H3)

The flow in $S_{\varepsilon}$ converges to the slow flow, as $\varepsilon\to 0$ .
(H4)

The set $S_{\varepsilon}$ is $C^{r}$ –smooth.
(H5)

The set $S_{\varepsilon}$ is normally hyperbolic and has the same stability properties with respect to the fast variables as $S_{0}$ .
(H6)

The set $S_{\varepsilon}$ is usually not unique. In regions that remain at a fixed distance from the topological boundary $\partial S_{\varepsilon}$ of $S_{\varepsilon}$ , all manifolds satisfying (H1)–(H5) lie at a Hausdorff distance $\mathcal{O}\left(e^{-K/\varepsilon}\right)$ from each other for some $K>0,$ where $K=\mathcal{O}\left(1\right)$ .

The proof of this theorem is extensive and covers topics in perturbation theory that are beyond the scope of this work. The reader may find it in [24] or in its original version in [7].

Definition 4.1.

The set $S_{\varepsilon}$ mentioned in Theorem 1 is referred to as the slow manifold of (3).

Discussion of Fenichel’s theorem

Theorem 1 is central to the study of systems with multiple time scales. Let us examine the role of each statement separately.

Statement (H6) is useful to simplify the language. In fact, there is not just one slow manifold $S_{\varepsilon}$ . However, since all possible slow manifolds $S_{\varepsilon}^{j},\,j\in\mathbb{N},$ are at a distance $\mathcal{O}\left(e^{-K/\varepsilon}\right)\overset{\varepsilon\to 0}{\longrightarrow}0$ from each other, we can simplify the discussion and consider only one slow manifold $S_{\varepsilon}$ .

The existence of a one-to-one correspondence between $S_{\varepsilon}$ and $S_{0}\subset\mathcal{C}_{0}$ , (H1), is important for the other statements in the theorem. The set $S_{\varepsilon}$ , being locally invariant, means that trajectories do not leave the slow manifold except at its boundary. Statement (H2) guarantees that the smaller $\varepsilon$ is, the closer $S_{\varepsilon}$ is to $S_{0}$ . Thus, since $S_{\varepsilon}$ is locally invariant, the trajectories on $S_{\varepsilon}$ approach $S_{0}$ by (H3). Finally, the fact that $S_{\varepsilon}$ is continuously differentiable to the same order as $S_{0}$ , (H4), and is diffeomorphic to $S_{0}$ , (H1), allows us to assert that the two sets $S_{0}$ and $S_{\varepsilon}$ share the same stability. The set $S_{\varepsilon}$ being normally hyperbolic, (H5), ensures that the fast variables have the same stability properties of $S_{0}$ .

The manifold $S_{\varepsilon}$

In this section, we give a precise way to describe $S_{\varepsilon}$ . Consider the (1,1)-fast–slow system in the slow time scale:

(15)			$\displaystyle\varepsilon\,\dfrac{\mathrm{d}x}{\mathrm{d}\tau}=\varepsilon\dot{x}=f(x,y,\varepsilon),$
(15)			$\displaystyle\phantom{\varepsilon\,}\dfrac{\mathrm{d}y}{\mathrm{d}\tau}=\phantom{\varepsilon}\dot{y}=g(x,y,\varepsilon),$

with $0<\varepsilon\ll 1$ .

Let $S_{0}$ be a compact, normally attracting subset of $\mathcal{C}_{0}$ without equilibria ( $g\left(x,y,\varepsilon\right)\neq 0$ ). By the implicit function theorem, there exists a function $x=x_{0}\left(y\right)$ , for $y\in\left]a,b\right[$ , such that we define $S_{0}$ as:

(16)		$\displaystyle S_{0}$	$\displaystyle=\left\{\left(x,y\right)\in{{\mathbb{R}}^{2}}\mathrel{\mathop{\mathchar 58\relax}}f\left(x,y,0\right)=0\wedge y\in\left]a,b\right[\right\}$
(16)			$\displaystyle=\left\{\left(x,y\right)\in{{\mathbb{R}}^{2}}\mathrel{\mathop{\mathchar 58\relax}}x=x_{0}\left(y\right)\wedge y\in\left]a,b\right[\right\}.$

Theorem 2 (Fenichel, 1979 [7]. See also Theorem 11.1.1 in Kuehn, 2015 [24].).

Let $S_{0}$ be a compact normally hyperbolic submanifold of $\mathcal{C}_{0}$ . Then there exists a slow manifold $S_{\varepsilon}$ that is $\mathcal{O}\left(\varepsilon\right)$ -close to $S_{0}$ for $\varepsilon>0$ sufficiently small. Locally, $S_{\varepsilon}$ is represented as a graph of a smooth function $h(y,\varepsilon)$ :

S_{\varepsilon}=\left\{\left(x,y\right)\in{\mathbb{R}}^{2}\mathrel{\mathop{\mathchar 58\relax}}x=h(y,\varepsilon)\right\},

where the map $y\mapsto h(y,\varepsilon)\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}\to{\mathbb{R}}$ has the asymptotic expansion

h(y,\varepsilon)=h_{0}\left(y\right)+h_{1}\left(y\right)\varepsilon+h_{2}\left(y\right)\varepsilon^{2}+{\mathcal{O}}(\varepsilon^{3}).

From the previous result, we may conclude that $h_{0}(y)=x_{0}(y)$ . The curve $h_{0}\left(y\right)$ defines locally the set $S_{0}$ of the critical manifold $\mathcal{C}_{0}$ .

In what follows, our aim is to show how to compute the term $h_{1}(y)$ in the expansion. Taking into account that $g\left(x,y,\varepsilon\right)\neq 0$ , we write (15) as

(17)

\varepsilon\,\dfrac{\mathrm{d}x}{\mathrm{d}y}=\dfrac{f\left(x,y,\varepsilon\right)}{g\left(x,y,\varepsilon\right)}.

By invariance, the set $x\equiv h(y,\varepsilon)$ satisfies (17) and hence

(18)

\varepsilon\dfrac{dh_{0}}{dy}(y)+\varepsilon^{2}\dfrac{dh_{1}}{dy}(y)+{\mathcal{O}}(\varepsilon^{3})=\dfrac{f\left(h(y,\varepsilon),y,\varepsilon\right)}{g\left(h(y,\varepsilon),y,\varepsilon\right)}.

We want to estimate the first terms in (18). For the first one, if $\left(x\left(\tau\right),y\left(\tau\right)\right)$ is a solution of (15) for $\varepsilon=0$ with initial condition in $S_{0}$ , hence $f\left(x(\tau),y(\tau),0\right)\equiv 0$ , then

f_{x}\left(x,y,0\right)\dfrac{dx}{d\tau}+f_{y}\left(x,y,0\right)\dfrac{dy}{d\tau}=0\quad\Rightarrow\quad\dfrac{dx}{dy}=-\dfrac{f_{y}\left(x,y,0\right)}{f_{x}\left(x,y,0\right)}

which implies (because $x\equiv h(y,\varepsilon)$ )

(19)

\dfrac{\partial h}{\partial y}(y,0)=\dfrac{dh_{0}}{dy}(y)=-\dfrac{f_{y}\left(h_{0}(y),y,0\right)}{f_{x}\left(h_{0}(y),y,0\right)}.

Now we compute $h_{1}(y)$ . To do this, let $H(x,y,\varepsilon)\mathrel{\mathop{\mathchar 58\relax}}=\dfrac{f(x,y,\varepsilon)}{g(x,y,\varepsilon)}$ , so the right hand side of (18) is

(20)

H\left(h(y,\varepsilon),y,\varepsilon\right)=\varepsilon\left(\dfrac{\partial H}{\partial x}\left(h(y,0),y,0\right)h_{1}(y)+\dfrac{\partial H}{\partial\varepsilon}\left(h(y,0),y,0\right)\right)+{\mathcal{O}}(\varepsilon^{2}),

since $H\left(h(y,0),y,0\right)=0$ and $\displaystyle\frac{\partial h}{\partial\varepsilon}(y,0)=h_{1}(y)$ .

Next, we use again $f(h(y,0),y,0)=0$ to compute the derivatives of $H$ appearing in (20):

\dfrac{\partial H}{\partial x}\left(h(y,0),y,0\right)=\dfrac{f_{x}\left(h(y,0),y,0\right)}{g\left(h(y,0),y,0\right)},\qquad\dfrac{\partial H}{\partial\varepsilon}\left(h(y,0),y,0\right)=\dfrac{f_{\varepsilon}\left(h(y,0),y,0\right)}{g\left(h(y,0),y,0\right)}\\

and (20) becomes

(21)

H\left(h(y,\varepsilon),y,\varepsilon\right)=\varepsilon\left(\dfrac{f_{x}\left(h_{0}(y),y,0\right)}{g\left(h_{0}(y),y,0\right)}h_{1}(y)+\dfrac{f_{\varepsilon}\left(h_{0}(y),y,0\right)}{g\left(h_{0}(y),y,0\right)}\right)+{\mathcal{O}}(\varepsilon^{2}).

Substituting (19) and (21) into (18) yields

-\varepsilon\dfrac{f_{y}\left(h_{0}(y),y,0\right)}{f_{x}\left(h_{0}(y),y,0\right)}+{\mathcal{O}}(\varepsilon^{2})=\varepsilon\left(\dfrac{f_{x}\left(h_{0}(y),y,0\right)}{g\left(h_{0}(y),y,0\right)}h_{1}(y)+\dfrac{f_{\varepsilon}\left(h_{0}(y),y,0\right)}{g\left(h_{0}(y),y,0\right)}\right)+{\mathcal{O}}(\varepsilon^{2})

hence

(22)

h_{1}(y)=-\dfrac{f_{y}\left(h_{0}(y),y,0\right)g(h_{0}(y),y,0)}{\left(f_{x}\left(h_{0}(y),y,0\right)\right)^{2}}-\dfrac{f_{\varepsilon}\left(h_{0}(y),y,0\right)}{f_{x}\left(h_{0}(y),y,0\right)}.

Example 4.1.

Consider the system FH-N (LABEL:fhn). If we take $S_{0}\subset\mathcal{C}_{0\text{L}}$ , we get that $S_{0}$ is defined as

h_{0}\left(y\right)=-\dfrac{4\left(\dfrac{2}{3}\right)^{1/3}}{\left(9y+\sqrt{3\left(27y^{2}-256\right)}\right)^{1/3}}-\dfrac{\left(9y+\sqrt{3\left(27y^{2}-256\right)}\right)^{1/3}}{18^{1/3}}.

Additionally, using (22), we get

h_{1}\left(y\right)=\dfrac{h_{0}\left(y\right)-by-c}{\left(4-h_{0}^{2}\left(y\right)\right)^{2}}

and $S_{\varepsilon}$ is approximated by the graph of

h_{\varepsilon}\left(y\right)=h_{0}\left(y\right)\,+\,\dfrac{h_{0}\left(y\right)-by-c}{\left(4-h_{0}^{2}\left(y\right)\right)^{2}}\,\varepsilon.\qquad\diamond

Since we defined the slow manifold $S_{\varepsilon}$ as an asymptotic expansion starting at the critical manifold $\mathcal{C}_{0}$ , we can apply the same reasoning to determine the period of a periodic solution, $\mathcal{T}_{\gamma_{\varepsilon}}$ , of a fast-slow system, with $0<\varepsilon\ll 1$ (cf. Figure 4(b)). In particular, using the same argument used to prove (13) we have:

\mathcal{T}_{\gamma_{\varepsilon}}=\mathcal{T}_{\gamma_{0}}+\mathcal{O}\left(\varepsilon\right).

5. Bifurcations

To understand the effect of the parameters $b$ and $c$ , we analyse two cases of (LABEL:fhn), where the first has $b=0$ and $c\in{\mathbb{R}}$ , and the second has $b\in{\mathbb{R}}\setminus{\left\{0\right\}}$ and $c=0$ .

(i) For $c\in{\mathbb{R}}$ and $b=0$ , system (LABEL:fhn) may be recast into the form:

(23)			$\displaystyle\dfrac{\mathrm{d}x}{\mathrm{d}t}=x^{\prime}=-y+4x-x^{3}$
(23)			$\displaystyle\dfrac{\mathrm{d}y}{\mathrm{d}t}=y^{\prime}=\varepsilon(x-c),$

with $0<\varepsilon\ll 1$ . It has only one equilibrium $\mathrm{E}=\left(c,\varphi(c)\right)$ for $\varphi(x)\coloneqq 4x-x^{3}$ , which results from the intersection of the cubic curve $y=\varphi(x)$ and the vertical line $x=c$ .

In order to evaluate the stability of $\mathrm{E}$ , we compute the jacobian matrix $J_{E}$ of the vector field associated to (23) at $E$ and find its eigenvalues, say $\lambda_{\mathrm{E}}^{\pm}\left(c,\varepsilon\right)$ . Thus, we get

\displaystyle J_{\mathrm{E}}=\begin{bmatrix}\,4-3c^{2}\,\phantom{aaa}\,&-1\\ \varepsilon&0\end{bmatrix}\qquad\text{and}\qquad\lambda_{\mathrm{E}}^{\pm}\left(c,\varepsilon\right)=\dfrac{\left(4-3c^{2}\right)\pm\sqrt{\left(4-3c^{2}\right)^{2}-4\varepsilon}}{2}.

It is straightforward to conclude that:

(1)

The equilibrium $\mathrm{E}$ is an unstable node/focus if $|c|<\dfrac{2}{\sqrt{3}}$ and
(2)

$\mathrm{E}$ is a stable node/focus if $|c|>\dfrac{2}{\sqrt{3}}$ .

When $c=\pm\mathrm{c^{H}}\text{, with }\mathrm{c^{H}}=2/\sqrt{3}$ , and $0<\varepsilon\ll 1$ , we get $\lambda_{\mathrm{E}}^{\pm}\left(\mathrm{c},\varepsilon\right)=\pm\mathrm{i}\sqrt{\varepsilon}$ . For $\varepsilon>0$ , as shown in Figure 5, the stability of the equilibrium point in the neighbourhood of $\mathrm{c^{H}}$ changes from a stable focus to an unstable focus as $c$ decreases. Similarly, by symmetry, the dual stability transition occurs in the neighbourhood of $-\mathrm{c^{H}}$ as $c$ increases.

For $c=1.15<\mathrm{c^{H}}$ and $\varepsilon=1$ , the equilibrium point $\mathrm{E}$ is an unstable focus, and thus it is expected that trajectories starting close to the fold point $P_{+}$ diverge, for $t>0$ . Numerically, we found a periodic solution when $c\overset{<}{\approx}\mathrm{c^{H}}$ (cf. Figure 6(d)). This dynamics is typical of a subcritical Hopf bifurcation, with a stable periodic solution as described in the next result. Following the terminology of [12, Cap. IV, Section 2], subcritical means that the periodic solution exists for $c<c^{\mathrm{H}}$ . If the periodic solutions are found for $c>c^{\mathrm{H}}$ , the bifurcation would be called supercritical.

Theorem 3 (Hopf, 1942. See also Marsden & McCraken, 1976 [26] for English translation in Section 5).

Given a one-parameter family of differential equations in ${\mathbb{R}}^{2}$ with parameter $c$ :

	$\displaystyle\dot{x}$	$\displaystyle=f(x,y;c)$
	$\displaystyle\dot{y}$	$\displaystyle=g(x,y;c),\qquad f,g\in{{C^{k}}}\quad\mathrm{(k}\geq 2\mathrm{)},$

such that there exists an equilibrium point $\left(x_{0}\left(c\right),y_{0}\left(c\right)\right)\in{{\mathbb{R}}^{2}}$ for all $c\in{\mathbb{R}}$ and the eigenvalues of the jacobian matrix $J_{\left(x_{0}\left(c\right),y_{0}\left(c\right)\right)}$ can be written as $\lambda\left(c\right)=\alpha\left(c\right)\pm i\,\beta\left(c\right)$ , if for some $c^{*}\in{\mathbb{R}}$ the following holds:

\alpha\left(c^{*}\right)=0,\quad\beta\left(c^{*}\right)\neq 0\quad\text{ and }\quad{\dfrac{\mathrm{d}\alpha}{\mathrm{d}c}}({c^{*}})\neq 0,

then, for $c$ close to $c^{*}$ , there exists a periodic solution (limit cycle) around $\left(x_{0}\left(c\right),y_{0}\left(c\right)\right)$ .

System (23) satisfies the conditions of Theorem 3, thus confirming the numerical finding of Figure 6(d). The direction of bifurcation and the stability of the periodic solution will be discussed in more detail in Section 6 below.

(ii) For $b\in{\mathbb{R}}\setminus{\left\{0\right\}}$ and $c=0$ , system (LABEL:fhn) may be written as:

(24)			$\displaystyle\dfrac{\mathrm{d}x}{\mathrm{d}t}=x^{\prime}=-y+4x-x^{3}$
(24)			$\displaystyle\dfrac{\mathrm{d}y}{\mathrm{d}t}=y^{\prime}=\varepsilon(x-by),$

with $0<\varepsilon\ll 1$ . For all $b\in{\mathbb{R}}\setminus{\left\{0\right\}}$ , the origin, $\mathrm{E}_{0}=\left(0,0\right)$ , is an equilibrium point the system (24). If either $b<0$ or $b>1/4$ , the system has two symmetric equilibria at $x_{\pm}=\pm\sqrt{4-b^{-1}}$ , as a result of a Pitchfork bifurcation, as the reader may check in Figure 8(a). They are:

\mathrm{E}_{\pm}=\left(x_{\pm},\varphi(x_{\pm})\right),\quad\text{with}\quad\varphi(x)\coloneqq 4x-x^{3}.

For $b<0$ , the equilibrium $E_{0}$ is an unstable node/focus and $E_{-}$ and $E_{+}$ are saddles. For $0<b\leq 1/4$ , there is a stable limit cycle (similar scenario to that of (v) represented in Figure 3(a)) and only the origin is an equilibrium point of the system, but it still exhibits the same unstable node/focus dynamics.

The equilibrium $\mathrm{E}_{0}$ becomes a saddle when $b>1/4$ . Moreover, at $b=1/4$ two new equilibria $\mathrm{E}_{-}$ and $\mathrm{E}_{+}$ are created (cf. Figure 8(a)). Given that $E_{+}$ and $E_{-}$ are symmetric, the stability of both equilibrium points is the same. Thus, we will only study the stability of $\mathrm{E}_{+}$ by computing

\displaystyle J_{\mathrm{E}_{+}}=\begin{bmatrix}3\,b^{-1}-8\phantom{aaa}&-1\\ \varepsilon&-\varepsilon\,b\end{bmatrix},\quad\operatorname{tr}\left(J_{\mathrm{E}_{+}}\right)=-\varepsilon\,b+3\,b^{-1}-8\quad\text{and}\quad\det\left(J_{\mathrm{E}_{+}}\right)=2\varepsilon\left(4b-1\right),

where $\operatorname{tr}$ and $\det$ denote the trace and determinant operators. Note that $\det\left(J_{\mathrm{E}_{+}}\right)>0$ , since $b>1/4$ . The eigenvalues of $J_{\mathrm{E}_{+}}$ satisfy

\lambda_{\varepsilon}^{\pm}\left(b\right)=\dfrac{\operatorname{tr}\left(J_{\mathrm{E}_{+}}\right)\pm\sqrt{\operatorname{tr}^{\,2}\left(J_{\mathrm{E}_{+}}\right)-4\,\det\left(J_{\mathrm{E}_{+}}\right)}}{2}.

If we take

\mathrm{b}_{\varepsilon}^{\mathrm{H}}=\dfrac{-4+\sqrt{16+3\varepsilon}}{\varepsilon}>\dfrac{1}{4},\qquad 0<\varepsilon\ll 1

then at $b=\mathrm{b}_{\varepsilon}^{\mathrm{H}}$

\operatorname{tr}\left(J_{\mathrm{E}_{+}}\right)=0\quad\text{and}\quad\mathinner{\dfrac{\mathrm{d}\mathrm{Re}\left(\lambda_{\varepsilon}^{+}\right)}{\mathrm{d}\,b}\Biggr{\rvert}}_{b=\mathrm{b}_{\varepsilon}^{\mathrm{H}}}<0

and therefore, by Theorem 3, there is a Hopf bifurcation. In the limit, as $\varepsilon\to 0$ , we have $\mathrm{b}_{\varepsilon}^{\mathrm{H}}\to 3/8$ (cf. Figure 7).

Moreover, for $b=3/8$ , the fold points are equilibria of (24).

The numerical simulations shown in Figure 8(c) indicate that the Hopf bifurcation at $\mathrm{E}_{+}$ is supercritical, hence generating an unstable periodic solution. This will be confirmed analytically in Section 6 below. This means that all trajectories starting in the regions bounded by the periodic solutions are attracted to the equilibrium point $\mathrm{E}_{-}$ or $\mathrm{E}_{+}$ , also lying in the interior of that region. On the other hand, when trajectories start outside these regions, they converge to the global stable limit cycle which already existed when the equilibria were unstable. As observed in the example of Figure 8(c), the length (a numerical estimate of the perimeter of the curve) of the unstable periodic solutions grows very quickly for small variations in $b$ , but reaches an unexpected maximum for $b=b_{\varepsilon}^{Hom}$ . This happens when the periodic solution captures the saddle equilibrium $\mathrm{E}_{0}$ at a homoclinic bifurcation [15, 29].

Definition 5.1.

Consider a system of ordinary differential equations $\dfrac{\mathrm{d}X}{\mathrm{d}t}=F(X)$ , where $F\in C^{2}$ , $X\in{\mathbb{R}}^{n}$ for some $n\in\mathbb{N}$ . Assume that $X_{0}\in{\mathbb{R}}^{n}$ is an equilibrium and $\gamma$ is a non-stationary solution of the system. The trajectory/orbit $\gamma$ is said to be homoclinic to $X_{0}$ if

\gamma\left(t\right)\to X_{0}\quad\text{ as }\quad t\to\pm\infty.

The existence of a saddle point equilibrium is not a sufficient condition for the existence of a homoclinic trajectory, as this also requires the existence of an intersection between the stable and unstable manifolds of the equilibrium point.

We already knew that (24) has the saddle equilibrium $\mathrm{E}_{0}=(0,0)$ and we have observed numerically that, for $b=b_{\varepsilon}^{Hom}$ , there is a homoclinic trajectory to the origin $\mathrm{E}_{0}$ . In Figures 8(c) and 9, which consider the system (24) for $\varepsilon=0.5$ , the homoclinic orbit appears for $b=b_{\varepsilon}^{Hom}\approx 0.36932$ , very close to $b_{\varepsilon}^{H}\approx 0.3666$ . The bifurcation where a periodic solution is destroyed at a homoclinic orbit is called a homoclinic bifurcation [15, 29].

6. Singularities and Canards

As we have seen in case (i) of Section 5, when the system’s equilibrium is also a fold point, a subcritical Hopf bifurcation occurs, giving rise to a stable periodic solution (cf. Figure 6). In Sections 3 and 4, we found another periodic solution which is not related to the Hopf bifurcation. Thus, for different parameters, we have two stable periodic solutions with distinct dynamical properties. In this section, the canard phenomenon bridges the local periodic solution from the Hopf bifurcation to the global limit cycle.

As illustrated in Figure 6(d), the existence in (23) of a periodic solution around the equilibrium point implies that the orbit traverses a part of its trajectory close to the repelling region of the critical manifold $\mathcal{C}_{0}$ . This behaviour seems to be contradictory since we found a solution that is close to the unstable branch of $\mathcal{C}_{0}$ .

Definition 6.1.

A slow trajectory that is $\mathcal{O}\left(\varepsilon\right)$ away from the repelling region of the slow manifold for a time of order $\mathcal{O}\left(1\right)$ , is called a canard.

Canards are related to equilibria of the slow equation that occur at fold points of the slow manifold.

Definition 6.2.

Given a singular $(1,1)$ –fast-slow system parametrised by $\lambda\in{\mathbb{R}}$ . Let $\mathfrak{p}=\left(x,y\right)$ be a fold point. The point $\mathfrak{p}$ is called a singular fold if:

(25)			$\displaystyle f\left(\mathfrak{p},\lambda,0\right)=0,$		$\displaystyle f_{x}\left(\mathfrak{p},\lambda,0\right)=0,$
(25)			$\displaystyle f_{xx}\left(\mathfrak{p},\lambda,0\right)\neq 0,$		$\displaystyle f_{y}\left(\mathfrak{p},\lambda,0\right)\neq 0\quad\text{and}\quad g\left(\mathfrak{p},\lambda,0\right)=0.$

The definition of singular fold point differs from that of a fold point since it is required to coincide with an equilibrium point of system (3) (i.e., $g\left(\mathfrak{p},\lambda,0\right)=0$ ).

Definition 6.3.

A singular fold $\mathfrak{p}$ is said to be regular if:

(26)

g_{x}\left(\mathfrak{p},\lambda,0\right)\neq 0\quad\text{ and }\quad g_{\lambda}\left(\mathfrak{p},\lambda,0\right)\neq 0.

The behaviour around a regular singular fold is described in the next theorem, a concatenation of two important results, whose proof can be found in [21, 22, 24].

Theorem 4 (Krupa & Szmolyan, 2001a, 2001b [21, 22]. See also Theorems 8.1.3 and 8.2.1 of [24]).

Consider a (1,1)–fast-slow system where $\left(x,y\right)=\left(0,0\right)$ is a regular singular fold point for $\lambda=0$ , in the following normal form:

(27)			$\displaystyle x^{\prime}=-y\,l_{1}\left(x,y,\lambda,\varepsilon\right)+x^{2}\,l_{2}\left(x,y,\lambda,\varepsilon\right)+\varepsilon\,l_{3}\left(x,y,\lambda,\varepsilon\right)$
(27)			$\displaystyle y^{\prime}=\varepsilon\left(\pm x\,l_{4}\left(x,y,\lambda,\varepsilon\right)-\lambda\,l_{5}\left(x,y,\lambda,\varepsilon\right)+y\,l_{6}\left(x,y,\lambda,\varepsilon\right)\right).$

Assume that, for $\varepsilon=0$ , there is a slow trajectory connecting the repelling and attracting regions of the critical manifold $\mathcal{C}_{0}$ . Then, there exist $\varepsilon_{0}>0$ and $\lambda_{0}>0$ such that for $0<\varepsilon<\varepsilon_{0}$ and $|\lambda|<\lambda_{0}$ , the system has an equilibrium point $\mathfrak{p}\in{\mathbb{R}}^{2}$ near the origin where $\mathfrak{p}\to\left(0,0\right)$ as $\left(\lambda,\varepsilon\right)\to\left(0,0\right)$ .

Moreover, if $\mathfrak{p}$ is stable for $\lambda<0$ , there exists a continuous function $\lambda_{c}\mathrel{\mathop{\mathchar 58\relax}}\left[0,\varepsilon_{0}\right]\to{\mathbb{R}}$ that associates each value of $\varepsilon\in\left[0,\varepsilon_{0}\right]$ to a value $\lambda$ that gives rise to a canard in the system, asymptotically defined by:

\lambda_{c}\left(\varepsilon\right)=-\left(B+A\right)\varepsilon+\mathcal{O}\left(\varepsilon^{3/2}\right),

and there exists a continuous function $\lambda_{\mathrm{H}}\mathrel{\mathop{\mathchar 58\relax}}\left[0,\varepsilon_{0}\right]\to{\mathbb{R}}$ that associates each value of $\varepsilon\in\left[0,\varepsilon_{0}\right]$ to a value $\lambda$ that gives rise to a Hopf bifurcation in the system, asymptotically defined by:

\lambda_{\mathrm{H}}\left(\varepsilon\right)=-B\varepsilon+\mathcal{O}\left(\varepsilon^{3/2}\right),

where

A=\dfrac{-\dfrac{\partial l_{1}}{\partial x}+3\dfrac{\partial l_{2}}{\partial x}-2\dfrac{\partial l_{4}}{\partial x}+2l_{6}}{8}\qquad B=\dfrac{\dfrac{\partial l_{3}}{\partial x}+l_{6}}{2},

where the functions $l_{i},\,i=1,\ldots,6$ and their partial derivatives are evaluated at the point $\left(x,y,\lambda,\varepsilon\right)=\left(0,0,0,0\right)$ .

For $A\neq 0$ , the Hopf bifurcation is non-degenerate. The Hopf bifurcation is supercritical if $A<0$ , and subcritical if $A>0$ .

Whenever a system with a regular singular fold exists, there is a change of coordinates that transforms the system into (27). As noted in Section 4, the existence of canards is coupled with the existence of a Hopf bifurcation, which in such systems is referred to as a singular Hopf bifurcation. That is, the values of $\lambda_{c}$ that give rise to canards are at most $\mathcal{O}\left(\varepsilon\right)$ away from the values $\lambda_{\mathrm{H}}$ where a Hopf bifurcation occurs,

\lambda_{\mathrm{H}}-\lambda_{c}=A\,\varepsilon+\mathcal{O}\left(\varepsilon^{3/2}\right).

An important consideration is the empirical difficulty in finding canards, since the smaller $\varepsilon$ is, the narrower the interval of $\lambda$ values for which canards appear. This small interval separates the $\lambda$ values at which a limit cycle occurs from those where it either undergoes a Hopf bifurcation or does not have any periodic solution. Figure 6(a) and its magnified Figures 6(b) and 6(c) illustrate this phenomenon. For this reason, the dynamics produced in the system around this small interval is called the canard explosion.

Example 6.1.

Consider the case (i) of Section 5. Since at the Hopf bifurcation point $c^{\mathrm{H}}$ the equilibrium $\mathrm{E}=\left(c^{\mathrm{H}},\varphi\left(c^{\mathrm{H}}\right)\right)=\left(\dfrac{2}{\sqrt{3}},\dfrac{16}{3\sqrt{3}}\right)$ is a regular fold singularity, to determine the stability of the periodic solution, we need to transform the system into its normal form (see [15, 24] for more details), which involves a translation of the equilibrium and of the bifurcation parameter to the origin. More specifically, we consider the change of coordinates:

\left(\bar{x},\bar{y},\bar{c}\right)=\left(x-c^{\mathrm{H}},y-\varphi\left(c^{\mathrm{H}}\right),c^{\mathrm{H}}-c\right),

where $\varphi\left(c^{\mathrm{H}}\right)=4c^{\mathrm{H}}-\left(c^{\mathrm{H}}\right)^{3}$ . Omitting the bars $\bar{x}$ and $\bar{y}$ this gives rise to the following system equivalent to (23):

(28)			$\displaystyle x^{\prime}=-y+x^{2}\left(-2\sqrt{3}-x\right)$
(28)			$\displaystyle y^{\prime}=\varepsilon\left(x+\bar{c}\right).$

Since the parameters $c$ and $\bar{c}$ have opposite orientations, then the existence of a supercritical Hopf bifurcation with periodic solutions for $\bar{c}>0$ implies the existence of a subcritical Hopf bifurcation with periodic solutions for $c<c^{H}$ in the original system (23).

It is easy to verify that $\mathfrak{p}=\left(-\bar{c},\varphi\left(-\bar{c}\right)\right)\,\overset{\bar{c}\to 0}{\longrightarrow}\left(0,0\right)\,$ is an equilibrium of the system and that $\mathfrak{p}$ is stable for $\bar{c}<0$ . Additionally, we have

f_{x}\left(\mathfrak{p},0\right)=0,\quad f_{xx}\left(\mathfrak{p},0\right)=-4\sqrt{3}\neq 0\quad\text{and}\quad f_{y}\left(\mathfrak{p},0\right)=-1\neq 0,

and also,

g_{x}\left(\mathfrak{p},0\right)=1\neq 0\quad\text{and}\quad g_{\bar{c}}\left(\mathfrak{p},0\right)=1\neq 0,

so $\mathfrak{p}$ is a regular singular fold point for $\bar{c}=0$ . Note that the system is defined in the normal form required by Theorem 4 and we have:

l_{1}=1,\quad l_{2}=-2\sqrt{3}-x,\quad l_{3}=0,\quad l_{4}=1,\quad l_{5}=-1,\quad l_{6}=0\quad\text{and}\quad A=-3/8<0.

Thus, the Hopf bifurcation is supercritical for $\bar{c}$ (subcritical for $c$ ) and occurs when

\bar{\mathrm{c}}_{\mathrm{H}}\left(\varepsilon\right)=0\cdot\,\varepsilon+\mathcal{O}\left(\varepsilon^{3/2}\right)

and the canards occur for

\bar{\mathrm{c}}_{c}\left(\varepsilon\right)=-\dfrac{3}{8}\,\varepsilon+\mathcal{O}\left(\varepsilon^{3/2}\right).\qquad\diamond

Example 6.2.

Now consider the case (ii) of Section 5 (with $c=0$ ) and the equilibrium $E_{+}=(x_{+},x_{+}/b)$ , with $x_{+}=\sqrt{4-b^{-1}}$ . Then $E_{+}$ is a regular singular fold point for $b=3/8$ with $x_{+}=4/3$ and there is a Hopf bifurcation point at $b=b_{\varepsilon}^{H}$ with $b_{\varepsilon}^{H}\to 3/8$ as $\varepsilon\to 0$ . In order to apply Theorem 4 we change variables and parameters by

(\bar{x},\bar{y},\lambda)=\left(x-x_{+},y-8x_{+}/3,b-3/8\right)

to obtain the equivalent system

(29)

\begin{array}[]{lcl}\bar{x}^{\prime}&=&-\bar{y}-\bar{x}^{3}-3\bar{x}^{2}x_{+}\\ \bar{y}^{\prime}&=&\varepsilon\left(\bar{x}-\lambda\bar{y}-\frac{3}{8}\bar{y}\right).\end{array}

In the notation of Theorem 4 we have

l_{1}=-1\qquad l_{2}=-3x_{+}-\bar{x}\qquad l_{3}=0\qquad l_{4}=1\qquad l_{5}=1\qquad l_{6}=-\frac{3}{8}

hence $A=-\dfrac{1}{8}\left(3+\dfrac{3}{4}\right)<0$ therefore the bifurcation is supercritical for $\lambda$ , and also for $b$ in the original system (24). Since for $b<b_{\varepsilon}^{H}$ the eigenvalues of the matrix $J_{E_{+}}$ are positive, then the bifurcating periodic solution is unstable.

Figure 10 illustrates the emergence of canards in system (23) with $\varepsilon=0.5$ and $\varepsilon=0.1$ . As we may see, the stable periodic solutions emerge from the subcritical Hopf bifurcation for $c=2/\sqrt{3}$ (analogous for $c=-2/\sqrt{3}$ ). As mentioned above, some of these periodic solutions remain close to the repelling region of the critical manifold, making them canards. In Figures 10(a) and 10(c), two types of orbits can be distinguished: the first type, coloured in yellow, consists of slow trajectories near $\mathcal{C}_{0\text{M}}$ and $\mathcal{C}_{0\text{R}}$ and one fast trajectory; the second type, coloured in green, which surrounds the other fold point, consists of three slow trajectories near $\mathcal{C}_{0\text{L}}$ , $\mathcal{C}_{0\text{M}}$ and $\mathcal{C}_{0\text{R}}$ and two fast trajectories. We refer to these two types of canards as headless canard and headed canard, respectively²²2This terminology is motivated by the literal meaning of the word “canard” in French (canard = duck). In this system, it is a bit more complicated to understand why, but if the reader tilts their head to the left they will see that the canard orbits resemble the shape of a duck, with or without head..

Note that the length of the canards grows almost instantaneously. As stated in [24], this was expected, but it is still surprising. The nearly vertical lines observed in Figures 10(b) and 10(d) show canards with lengths, $\mathcal{A}\left(c\right)\in\left[2,20\right]$ for $c\approx 1.150077$ in the case $\varepsilon=0.5$ and for $c\approx 1.153794$ in the case $\varepsilon=0.1$ .

7. Discussion and Further Work

This article studies the dynamics of a fast-slow system derived from the FH-N model according to the geometric singular perturbation theory and the bifurcation theory methods. In particular, we provide an analytic proof that the fast-slow FH-N system presents relaxation oscillation dynamics as well as periodic solutions induced by Hopf bifurcation. We emphasise the emergence of a homoclinic orbit and canards connecting these two distinct dynamics. We illustrate each result with numerical computations. In each section, we have complemented the discussion of observed dynamics with stability and bifurcation analysis. This approach empowers the readers to navigate the parameter space, enabling them to find specific dynamical regimes of interest.

The FH-N equation is an excellent candidate to serve as a starting point for developing new methods to analyse transient and non-reciprocal coupled interactions. Cardiac and neuronal examples, in particular, highlight the FH-N model’s adaptability and significance in these areas [2].

When coupling two FH-N systems through the slow equations, we have performed numerical simulations that suggest the existence of mixed-mode oscillations where the dynamics is characterised by the existence of orbits that alternate periodically between large and small amplitude oscillations [13]. Recently, [20] showed that the coupling through the fast equations of two distinct FH-N models produces mixed-mode oscillations induced by a cusp singularity present in this system. We believe that a similar result may explain our numerical findings.

The prevalence of chaotic attractors in periodically perturbed fast-slow systems is yet to be explored. The article [16] bridged the areas of geometric singular perturbation theory and of chaotic attractors theory and provided a general technique for proving the existence of chaotic attractors for periodically perturbed two-dimensional vector fields with two time scales, that are equivalent to:

(30)			$\displaystyle\varepsilon\,\dfrac{\mathrm{d}x}{\mathrm{d}t}=f(x,y,\theta)$
			$\displaystyle\phantom{\varepsilon\,}\dfrac{\mathrm{d}y}{\mathrm{d}t}=g(x,y,\theta)$
			$\displaystyle\phantom{\varepsilon\,}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}=\omega,$

where $(x,y,\theta)\in{\mathbb{R}}\times{\mathbb{R}}\times{\mathbb{S}}^{1}$ and $\omega>0$ is the slow driving frequency.

Results in [16] connect two areas of dynamical systems: the theory of chaotic attractors for discrete two-dimensional Hénon-like maps and geometric singular perturbation theory. Two-dimensional Hénon-like maps are diffeomorphisms whose singular limit is a family of non-invertible one-dimensional maps. In [32], the authors obtained sufficient conditions for the existence of chaotic attractors in these families. Three-dimensional singularly perturbed vector fields have return maps that are also two-dimensional diffeomorphisms with folds. To represent fully the behaviour of system (30) in the singular limit, we must allow jumps of trajectories from one sheet of the critical manifold to another that follow the direction of trajectories when $\varepsilon>0$ . Jumps parallel to the $x$ -axis occur at folds, where the tangent plane to the critical manifold includes this direction. We have chosen to leave the study of periodically perturbed FH-N models for future work.

Acknowledgments

The first and second authors were partially supported by CMUP, member of LASI, which is financed by national funds through FCT – Fundação para a Ciência e a Tecnologia, I.P., under the projects with reference UIDB/00144/2020 and UIDP/00144/2020. The third author has been supported by the Project CEMAPRE/REM – UIDB/05069/2020 financed by FCT/MCTES through national funds.

References

[1] M.-S. Abdelouahab, R. Lozi and G. Chen, Complex canard explosion in a fractional-order FitzHugh–Nagumo model, International Journal of Bifurcation and Chaos 29 1950111 (2019).
[2] D. Cebrián-Lacasa, P. Parra-Rivas, D. Ruiz-Reynés, and L. Gelens, Six decades of the FitzHugh-Nagumo model: A guide through its spatio-temporal dynamics and influence across disciplines (2024).
[3] S.–N. Chow, C. Li and D. Wang, Normal Forms and Bifurcation of Planar Vector Fields, Cambridge University Press (1994).
[4] M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H.M. Osinga and M. Wechselberger, Mixed-mode oscillations with multiple time scales, SIAM Review 54 211–288 (2012)
[5] X. Dong and C. Wang, Identification of the FitzHugh–Nagumo model dynamics via deterministic learning, International Journal of Bifurcation and Chaos 25, 1550159 (2015).
[6] C. Doss-Bachelet, J.-P. Françoise, and C. Piquet, Bursting oscillations in two coupled FitzHugh-Nagumo systems, ComPlexUs 1 (2003), no. 3, 101–111.
[7] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, Journal of Differential Equations 31 (1979), no. 1, 53–98.
[8] R. FitzHugh, Thresholds and plateaus in the Hodgkin–Huxley nerve equations, Journal of General Physiology 43 , no. 5, 867–896 (1960).
[9] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophysical journal 1 (1961), no. 6, 445–466.
[10] A. Georgescu, C. Rocşoreanu and N. Giurgiţeanu, Global Bifurcations in FitzHugh-Nagumo Model in Bifurcation, symmetry and patterns, Birkhäuser Basel, Basel 197–202 (2003).
[11] J.-M. Ginoux and J. Llibre, Canards existence in FitzHugh-Nagumo and Hodgkin-Huxley neuronal models, Mathematical Problems in Engineering (2015)
[12] M. Golubitsky and D.G. Schaeffer, Singularities and groups in bifurcation theory, Vol. 1, Springer–Verlag (1985).
[13] B.F.F. Gonçalves, I.S. Labouriau, and A.A.P. Rodrigues, Bifurcation and canards in coupled FitzHugh-Nagumo equations, In preparation, 2024.
[14] W. Govaerts, Yu. A. Kuznetsov, H. G. E. Meijer, B. Al-Hdaibat, V. De Witte, A. Dhooge, W. Mestrom, N. Neirynck, A. M. Riet and B. Sautois, Matcont: Continuation toolbox for odes in matlab, 2019.
[15] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, Springer–Verlag, 1983.
[16] J. Guckenheimer, M. Wechselberger and L.-S. Young, Chaotic attractors of relaxation oscillators, Nonlinearity 19 (2006), no. 3, 701.
[17] P. Hartman, Ordinary differential equations, 2nd ed., Classics in applied mathematics 38, Society for Industrial and Applied Mathematics, 1982.
[18] A. L. Hodgkin and A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, Journal of Physiology 117 (1952), no. 4, 500–544.
[19] J. Keener and J. Sneyd, Mathematical physiology, Springer–Verlag, 1998.
[20] K. U. Kristiansen and M. G. Pedersen, Mixed-mode oscillations in coupled FitzHugh-Nagumo oscillators: Blow-up analysis of cusped singularities, SIAM Journal on Applied Dynamical Systems 22 (2023), no. 2, 1383–1422.
[21] M. Krupa and P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions, SIAM Journal on Mathematical Analysis 33 (2001), no. 2, 286–314.
[22] M. Krupa and P. Szmolyan, Relaxation oscillation and canard explosion, Journal of Differential Equations 174 (2001), no. 2, 312–368.
[23] M. Krupa, M., A. Vidal, M. Desroches and F. Clément, Mixed-mode oscillations in a multiple time scale phantom bursting system, SIAM Journal on Applied Dynamical Systems 11, 1458–1498 (2012).
[24] C. Kuehn, Multiple time scale dynamics, Springer–Verlag, 2015.l
[25] J. Llibre and C. Vidal, Periodic solutions of a periodic FitzHugh-Nagumo system International Journal of Bifurcation and Chaos 25 (2015).
[26] J. E. Marsden and M. McCracken, The Hopf bifurcation and its applications, Springer-Verlag, 1976.
[27] J. Nagumo, S. Arimoto, and S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proceedings of the IRE 50 (1962), no. 10, 2061–2070.
[28] P.E. Phillipson and P. Schuster, A comparative study of the Hodgkin–Huxley and FitzHugh–Nagumo models of neuron pulse propagation, International Journal of Bifurcation and Chaos 15, 3851–3866 (2005).
[29] S. H. Strogatz, Nonlinear dynamics and chaos, CRC Press, 2015.
[30] S. Tsuji, T. Ueta, H. Kawakami, and K. Aihara, A design method of bursting using two-parameter bifurcation diagrams in FitzHugh–Nagumo model, International Journal of Bifurcation and Chaos 14 2241–2252 (2004).
[31] B. Wang, Dynamical behavior of the almost-periodic discrete FitzHugh–Nagumo systems, International Journal of Bifurcation and Chaos 17 1673–1685 (2007).
[32] Q. Wang and L.-S. Young, Toward a theory of rank one attractors, Annals of Mathematics 167 (2008), no. 2, 349–480.

Bifurcations and canards in the FitzHugh-Nagumo system: a tutorial in fast-slow dynamics