# Mixed-Mode Oscillations Based on Complex Canard Explosion in a Fractional-Order Fitzhugh-Nagumo Model

René Lozi, Mohammed-Salah Abdelouahab, Guanrong Chen

<figure><img src="image_1.png" /></figure>

## To cite this version:

René Lozi, Mohammed-Salah Abdelouahab, Guanrong Chen. Mixed-Mode Oscillations Based on Complex Canard Explosion in a Fractional-Order Fitzhugh-Nagumo Model. Applied Mathematics and Nonlinear Sciences, 2020, 5 (2), pp.239-256. 10.2478/AMNS.2020.2.00047. hal-02541322

HAL Id: hal-02541322

https://hal.science/hal-02541322

Submitted on 13 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mixed-Mode Oscillations Based on Complex Canard Explosion in a Fractional-Order Fitzhugh-Nagumo Model.

René Lozi†, Mohammed-Salah Abdelouahab, Guanrong Chen.

Laboratoire de Mathématiques J.A. Dieudonné, Université Côte d'Azur, Parc Valrose, 06108 Nice Cedex 02
France

Laboratory of Mathematics and Their Interactions, University Abdelhafid Boussouf, Mila
Algeria

Department of Electrical Engineering, City University of Hong Kong, Kowloon, Hong Kong SAR, P.R.
China

**Abstract**

This article highlights particular mixed-mode oscillations (MMO) based on canard explosion observed in a fractional-order Fitzhugh-Nagumo (FFHN) model. In order to rigorously analyze the dynamics of the FFHN model, a recently introduced mathematical notion, the Hopf-like bifurcation (HLB), which provides a precise definition for the change between a fixed point and an *S*-asymptotically *T*-periodic solution, is used. The existence of HLB in this FFHN model is proved and the appearance of MMO based on canard explosion in the neighborhoods of such HLB points are numerically investigated using a new algorithm: the global-local canard explosion search algorithm. This MMO is constituted of various patterns of solutions with an increasing number of small-amplitude oscillations when two key parameters of the FFHN model are varied simultaneously. On the basis of such numerical experiment, it is conjectured that chaos could occur in a two-dimensional fractional-order autonomous dynamical system, with the fractional-order close to one. Therefore, this very simple two-dimensional FFHN model, presents an incredible ability to mimic the complex dynamics of neurons.

**Keywords :** Fitzhugh-Nagumo model ; canard explosion ; mixed-mode oscillation.

**AMS 2010 codes :** 34A08, 34C15,34E17,37G15.

<footnote>†Corresponding author.</footnote>

<footnote>Email address : René.lozi@univ-cotedazur.fr</footnote>
# 1 Introduction

The most popular physical model (an electric circuit) mimicking the action potential of many types of neurons was introduced in 1952, by Hodgkin and Huxley [1] who were awarded the Nobel Prize in Physiology or Medicine in 1963 for this work.

In biophysics science, this model (thereafter called HH) is often regarded as one of the major achievements of the last-century. Since then, it has been improved in many ways, for example in [2] to incorporate transition state theory and produce thermodynamic HH model.

The Kirchhoff laws applied to this circuit afford a mathematical system of four Ordinary Differential Equations (ODE) controlled by at least 13 parameters. Although nowadays it is very easy to compute numerically the solutions of this system, due to the huge number of parameters and the four dimensions of the phase space, it is very difficult to study it thoroughly.

In order to cope with the challenging task of mimicking the action potential of many types of neurons, several other channel-based models comporting a great number of nonlinear equations have been designed to capture the physiological processes in the membrane (see [3, 4] for a survey).

Contrarily, phenomenological models, which try only to replicate the characteristic features of the bursting behavior without direct relation to what happens in the neuron at biological level, have been derived from HH model.

Among them, the most known are the Fitzhugh-Nagumo (FHN) and the Hindmarsh-Rose models. The FHN model [5] is a system of two-equation reduction of the system of four-equation system HH. There are two different Hindmarsh-Rose models: the first one was published in a paper dated 1982 [6], which comes from a modification of the FHN model, the second one [7] is a more sophisticated 3-D model, which is known to demonstrate almost all types of robust activities of the HH model.

In this article, we restrain ourselves to the FHN model. However, we consider its fractional version, because on one hand, the four-dimensional nonlinear differential equation of the HH model is difficult to study in details; on the other hand, due to the innermost properties of two-dimensional dynamical systems highlighted by the Poincaré-Bendixon theorem [8], the FHN model is unable to reproduce many complex dynamics of the corresponding four-dimensional system, such as chaos and hyperchaos.

Moreover, mixed-mode oscillations (MMO), which are very common in electroencephalography (EEG) data, can only be modeled by autonomous systems of ODE with integer dimensions greater than two. Nevertheless, it is possible to find a resolution to this concerned issue between too simple and too complex systems, using fractional derivatives.

Recently, some articles have partially studied the fractional-order Fitzhugh-Nagumo (FFHN) model [9–12], or its modified versions. In this article, we investigate the MMO based on complex canard explosion in the FFHN model. The nature of systems of fractional ODE prohibits the existence of periodic solutions on a finite time interval [13], therefore the existence of Hopf bifurcation. Consequently, to rigorously analyze the FFHN, a Hopf-like bifurcation (HLB) theory is applied [14], which provides a precise definition for the change between a fixed point and an *S*-asymptotically *T*-periodic solution of a fractional-order system. The emergence of MMO solutions formed by patterns of one large oscillation followed by several small-amplitude oscillations of the FFHN, when the fractional-order is close to one, is highlighted. In particular, when two parameters are varied simultaneously, one can observe that the number of such small-amplitude oscillations increases in every pattern. Such a phenomenon is forbidden from a system with the order of derivative being equal to one, due to the Poincaré-Bendixon theorem.

On can infer, from the study of the complex canard explosion and MMO in the FFHN model, that chaotic phenomenon can occur in two-dimensional autonomous fractional-order dynamical systems, with a fractional-order close to one. In conclusion, the FFHN model is a very simple two-dimensional model with an incredible ability to present the complex dynamics of neurons.

This article is organized as follows: in Sec. 2, the definition and properties of canard and MMO will be
reviewed. In Sec. 3, some classical results on fractional calculus will be recalled. In Sec. 4, fractional calculus will be further discussed within the context of dynamical systems. In Sec. 5, the new definition of HLB will be presented, followed by an analysis of the stability and HLB of the solutions of the FFHN model. In Sec. 6, complex canard explosion and MMO in the FFHN model will be investigated. Finally, in Sec. 7, a brief conclusion will be drawn.

## 2 Mixed-Mode Oscillations (MMO) Based on Canard solutions

On one hand, many electric or electronic devices, like oscilloscopes, display curves as results of sampled physical or physiological phenomena. In medicine, the shapes of these curves are often used to characterize diseases like epilepsy and stroke. On the other hand, biological models involving ODE provide curves in the phase space as solutions. Therefore, the study of particular patterns in such trajectory curves is very important in practice. Among these patterns, MMO combine features of small and large oscillations of relaxation type. It is in the well-known Belousov-Zhabotinsky chemical reaction that they were first discovered [15]. Some MMO are due to the existence of canard solutions [16]. Those solutions of slow-fast systems of differential equation are trajectories jumping from stable to unstable parts of a slow manifold, and vice versa.

### 2.1 Canard trajectories

Studying the van der Pol equation in the Liènard plane [17], a team of French mathematicians, habitually using the theory of nonstandard analysis (NSA), discovered a very strange bifurcation phenomenon [18]. They coined the French name of canard (duck) for such kind of bifurcation which occurs within an exponentially small range of parameters, highlighting the very stiff transition from a large-amplitude limit cycle (relaxation) to a small one in a slow-fast ODE, which is also referred to as singularly perturbed systems. Nowadays, the name canard, even not an English name, is routinely accepted by the mathematical community and the corresponding phenomenon is widely applied for various purposes without the use of NSA [19]. Its prototypical model is the van der Pol system with constant forcing $a > 0$ [20], described in the Liènard plane $(x, y)$ as

$$ \begin{cases} \dot{x} = y - \frac{1}{3}x^3 + x, \\ \dot{y} = \varepsilon(a - x). \end{cases} \qquad (1) $$

For a small positive parameter $\varepsilon \ll 1$, the variable $x$ is driven by the fast vector field evolving on the fast time scale $t = \tau/\varepsilon$, and $y$ evolves on a slow time scale $\tau$ directed by the slow vector field. Thus, $x$ (resp. $y$) is called the fast (resp. slow) variable. The cubic shaped curve $S$, defined by $\dot{x} = 0$, is called the nullcline, or slow curve, or the critical manifold. Among its three branches, the middle one is repelling and connected with both attracting outer branches via two fold-points, as shown in Fig. 1. Starting from any initial point in the Liènard plane, the solution of (1) is quickly attracted by the fast vector field in a neighborhood of one stable branch of the slow curve. As an example, if the initial point belongs to the upper right-hand side of this plane, the solution of (1) will follow this stable nullcline downward, directed by the slow vector field, until it reaches the lower fold-point from which it jumps leftward to the other stable part of the critical manifold, as shown in Fig. 1(a). The slow vector field will then drive the trajectory upward until it reaches the other fold-point. Then, the trajectory will jump rightward to the previous stable part of the critical manifold, forming a periodic cycle. However, it is possible that, depending on the value of the parameter $a$, instead of jumping immediately from the bottom fold-point to the other stable slow curve, the trajectory follow the unstable part for a short time, as shown in Figs. 1(b) and 1(c). This “amplitude bifurcation”, in which the periodic nature of the solution is unchanged, is referred to as canard, because its shape resembles that of a duck [Fig. 1(b), in which both legs are maliciously added].

In general, canard solutions occur in singularly perturbed systems of the form
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Canard explosion of the van der Pol oscillator for $\epsilon = 0.02$ happens in an exponentially small parameter interval near $a = a^* \approx 0.997461066999$ where the transition from relaxation oscillations for (a) $a < 0.997461066999$, to small-amplitude limit cycles for (c) $a \ge 0.997461066999$, comes via (b) canard cycles.</figcaption></figure>

$$ \left\{ \begin{aligned} \varepsilon \frac{dx(\tau)}{d\tau} &= f(x(\tau), y(\tau), \varepsilon), \\ \frac{dy(\tau)}{d\tau} &= g(x(\tau), y(\tau), \varepsilon), \end{aligned} \right. \qquad (2) $$

with $x \in \mathbb{R}^n$ and $y \in \mathbb{R}^m$ being the fast (resp. the slow) variable, both functions $f$ and $g$ are sufficiently smooth, and $\tau$ denotes the independent variable of (2), called the slow time scale. After a time re-scaling to the fast time scale $t = \tau/\epsilon$ (for $\epsilon \neq 0$), (2) is equivalent to the system

$$ \left\{ \begin{aligned} \dot{x} &= f(x, y, \varepsilon), \\ \dot{y} &= \varepsilon g(x, y, \varepsilon), \end{aligned} \right. \qquad (3) $$

where the dot denotes differentiation with respect to $t$.

As, $\epsilon \to 0$, system (2) is reduced to the slow subsystem

$$ \left\{ \begin{aligned} 0 &= f(x, y, 0), \\ \dot{y} &= g(x, y, 0), \end{aligned} \right. \qquad (4) $$

and system (3) to the fast (or layer) subsystem

$$ \left\{ \begin{aligned} \dot{x} &= f(x, y, 0), \\ \dot{y} &= 0. \end{aligned} \right. \qquad (5) $$

Both subsystems are not equivalent, but they play a key role in the geometric singular perturbation theory that deals with the dynamical analysis of the full system for $\epsilon > 0$ [21]. The critical (slow) manifold

$$ S = \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^m : f(x, y, 0) = 0\} $$

for system (2) corresponds to the phase space of the reduced subsystem (4) and the set of equilibria of the layer subsystem (5).

The subset $S^a \subset S$, for which all eigenvalues of $D_x f$ have negative real parts, is called the attracting slow invariant manifold and, similarly, the subset $S^r \subset S$ for which all eigenvalues of $D_x f$ have positive real parts is called the repelling slow invariant manifold. It is then possible to define rigorously a canard solution of (3).
**Definition 1.** [22] A solution of (2) or (3), first following $S^a$ and then continuing for a while along $S^r$, is called a canard solution. Conversely, the solution first follows $S^r$, and then continuing after $S^r$ for a while, is called a false canard solution.

## 2.2 Mixed-mode oscillations

MMO consists of $L$ large-amplitude (relaxation) oscillations followed by $s$ small-amplitude (subthreshold) oscillations, simply denoted by $L^s$ [16]. Canard phenomenon helps to understand and analyze such kind of slow fast dynamics. Complex oscillatory patterns like MMO can be explained by the coupling of local passage near a folded singularity, around which canard solutions emerge, with the global return mechanism via relaxation spikes that reset the local dynamics. The topic of MMO is of great importance in various applications, such as cellular electrical and secretory activities [23, 24], chemical reactions [25], optical oscillations [26], etc. In 1963, the Nobel Prize for Medicine was awarded to Alan Lloyd Hodgkin and Andrew Huxley for their work of 1952 when they built a mathematical model (Hodgkin-Huxley model) of electric circuits, which reproduces fairly accurately the action potential of many types of neurons. This model is a nonlinear system of four ODE. The complexity of the mathematical analysis of such system had motivated the introduction of various simplifications of it, the best known (and simplest) one being probably the Fitzhugh-Nagumo model proposed in 1961 [5], described by a two-dimensional differential system with cubic nonlinearity. Unfortunately, bounded solutions of a two-dimensional autonomous system are attracted to fixed points or period cycles. They cannot produce chaos and hyperchaos and all the complex dynamics of a four-dimensional dynamical system like the Hodgkin-Huxley model, proved by the Poincaré-Bendixon theorem [8]. Indeed, MMO dynamics can only be modeled by autonomous systems of ODE of greater than two dimensions. The occurrence of chaotic behavior in the Hodgkin-Huxley model was reported in [27] as well as the occurrence of the MMO in [28].

# 3 The Basis of Fractional Calculus

Although the roots of fractional differential calculus go back to the early days of differential calculus [29], still few mathematicians use it today, even if their number of publications increases from year to year. Fractional differential calculus can be understood in the scope of generalization of integration and differentiation. Many systems in interdisciplinary fields take an advantage to be described by fractional-order differential equations, such as viscoelastic systems, dielectric polarization, and quantum evolution of complex systems [30–33].

The most interesting property of fractional systems is that they allow to model complex phenomena, without having to increase their dimension as for classical systems. It is this property that we will use in Section 4, to model the Fitzhugh-Nagumo system [5], with only two equations.

Several definitions of fractional-order derivatives have been introduced in the literature [34–37]. A common one is the Riemann-Liouville definition of fractional derivatives [34, 35], given by

$$ 
\begin{aligned}
\frac{R}{a} D_t^\alpha x(t) &= \frac{1}{\Gamma(m-\alpha)} \frac{d^m}{dt^m} \int_a^t (t-\tau)^{m-\alpha-1} x(\tau) d\tau \\
&= \frac{d^m}{dt^m} (a j_t^{m-\alpha} x(t)), \quad t > a, m-1 \le \alpha < m,
\end{aligned}
$$

where $\Gamma$ is the gamma function and $a j_t^\beta$ is the Riemann-Liouville integral operator defined as

$$ a j_t^\beta x(t) = \frac{1}{\Gamma(\beta)} \int_a^t (t-\tau)^{\beta-1} x(\tau) d\tau. $$

The Laplace transform of the $\alpha$-order Riemann-Liouville differential operator is

$$ L\{\int_0^t D_t^\alpha x(t)\} = s^\alpha L\{x(t)\} - \sum_{k=0}^{m-1} s^k \int_0^t D_t^{\alpha-1-k} x(t) \Big|_{t=0}, $$
which involves initial conditions expressed in terms of fractional derivatives that have no clear physical interpretation. To avoid this problem, Caputo [37] introduced a new definition as follows:

$$ 
\begin{aligned}
c_a^\alpha D_t^\alpha x(t) &= \frac{1}{\Gamma(n-\alpha)} \int_a^t (t-\tau)^{n-\alpha-1} x^{(n)}(\tau) d\tau \\
&= j^{n-\alpha} \left( \frac{d^n}{dt^n} x(t) \right), \quad t > a,
\end{aligned}
$$

where $n = \lceil \alpha \rceil$ is the value of $\alpha$ rounded up to the nearest integer. The Laplace transform of the $\alpha$-order Caputo differential operator is

$$ L\{c_0^\alpha D_t^\alpha x(t)\} = s^\alpha L\{x(t)\} - \sum_{k=0}^{m-1} s^{\alpha-1-k} x^{(k)}(0). $$

Clearly it involves initial condition in terms of only integer derivatives. The Grünwald-Letnikov definition [34] is given by

$$ c_a^\alpha D_t^\alpha x(t) = \lim_{h \to 0} \frac{1}{h^\alpha} \sum_{k=0}^{k=\lfloor t \rfloor} (-1)^k \frac{\Gamma(\alpha+1)}{k! \Gamma(\alpha-k+1)} x(t-kh), \quad (6) $$

where $t > a$ and $\alpha$ is a positive real number. Its integral form is

$$ c_a^\alpha D_t^\alpha x(t) = \lim_{h \to 0} h^\alpha \sum_{k=0}^{k=\lfloor t \rfloor} \frac{\Gamma(\alpha+k)}{k! \Gamma(\alpha)} x(t-kh). \quad (7) $$

When $x$ is of class $C^m$, where $m-1 \le \alpha < m$, both Riemann-Liouville and Grünwald-Letnikov definitions are equivalent. Hence, to compute numerical approximations of the Riemann-Liouville fractional derivative the Grünwald-Letnikov definition is adopted.

# 4 Fractional Fitzhugh-Nagumo Model and Hopf-like Bifurcation

In this section, we focus on a fractional version of the FHN model in order to give it more flexibility, so that its solutions are closer to the initial model of Hodgkin-Huxley which it approximates, without increasing its number of equations.

## 4.1 The electrical circuit corresponding to the Fitzhugh-Nagumo model

An electrical circuit equivalent to the classical FHN model was constructed by Nagumo et al. [38], as shown in Fig. 2. It consists of a voltage variable $v$ (membrane potential) with cubic nonlinearity that allows regenerative self-excitation via a positive feedback, and a recovery variable $w$, which describes the combined effect of ion channels, with a linear term that affords a slower negative feedback.

The equation of this model reads

$$ \begin{cases} \frac{dv}{dt} = v - \frac{v^3}{3} - w + I \\ \frac{dw}{dt} = \frac{1}{T}(v + a - bw) \end{cases} \quad (8) $$

with parameters $I$, $a$, $b$ and $T$.

On can reduce it to the form of Eq. (1), by posing $x = v$, $y = w$ and $\varepsilon = 1/T$, as

$$ \begin{cases} \dot{x} = x - \frac{x^3}{3} - y + I, \\ \dot{y} = \varepsilon(x + a - by). \end{cases} \quad (9) $$

It is obvious that, if $b = I = 0$ and $x \to -x$, Eq. (9) is exactly the van der Pol system (1).
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Electrical circuit equivalent to the Fitzhugh–Nagumo model.</figcaption></figure>

## 4.2 The FHN model

In fact, the equivalence between the circuit in Fig. 2 and Eq. (9) is not exact, because, as shown by Jonscher [39], an ideal capacitor has an integer-order constitutive equation

$$I(t) = C \frac{dv(t)}{dt},$$

where $I(t)$ is the current through the capacitor and $v(t)$ the voltage across it, cannot exist in nature. In [40], it shows that a realistic capacitor could better be represented with a fractional-order constitutive equation,

$$I(t) = CD^{\alpha}v(t),$$

where $\alpha$ is a constant related to the proximity effect. Moreover, the first author had already showed [41] that an inductor is fractional in nature with the constitutive relationship between voltage and intensity:

$$v(t) = LD^{\alpha}I(t),$$

where the derivative order $\alpha$ is a constant related to the loss of the capacitor. Based on these considerations, Liu et al. [9] introduced the fractional-order version of the Fitzhugh–Nagumo model (FFHN):

$$\begin{cases} D^{\alpha_1}x = x - \frac{1}{3}x^3 - y + I, \\ D^{\alpha_2}y = \varepsilon(x + a - by), \end{cases} \quad (10)$$

where the derivative order $\alpha_1$ (resp. $\alpha_2$ is constant and related to the loss of the capacitor (resp. the proximity effect of the inductor).

## 4.3 Stability of the fixed point of FFHN

In this section, we briefly recall some results concerning the fixed points of the FFHN for $\alpha_1 = \alpha_2 = \alpha \in (0, 2)$. In order to have a unique fixed point, some parameter restrictions are introduced.
**Proposition 1.** [14] For all $a, b, I \in \mathbb{R}$ satisfying

$$-4 \left(1 - \frac{1}{b}\right)^3 + 9 \left(I - \frac{a}{b}\right)^2 > 0,$$

the system (10) has a unique equilibrium point

$$E = (x_e(b), y_e(b)),$$

where

$$x_e(b) = \sqrt[3]{\frac{-q + \sqrt{\frac{-\Delta}{27}}}{2}} + \sqrt[3]{\frac{-q - \sqrt{\frac{-\Delta}{27}}}{2}}, y_e(b) = x_e(b) - \frac{1}{3}x_e^3(b) + I, \\ \Delta = -(4p^3 + 27q^2), \quad p = -3\left(1 - \frac{1}{b}\right) \quad \text{and} \quad q = -3\left(I - \frac{a}{b}\right).$$

**Proposition 2.** [14] The fixed point E of Eq. (10) is locally asymptotically stable if

$$\left| \arctan \left( \frac{\sqrt{-(x_e^2(b) - b\varepsilon - 1)^2 + 4\varepsilon}}{x_e^2(b) + b\varepsilon - 1} \right) \right| > \alpha \frac{\pi}{2}. \quad (11)$$

The proof is straightforward based on the following theorems.

**Theorem 3.** [49, 50] The following fractional-order linear autonomous system:

$$\begin{cases} D^{\alpha} X = AX, \\ X(0) = X_0, \end{cases} \quad X \in \mathbb{R}^n, 0 < \alpha < 2 \text{ and } A \in \mathbb{R}^n \times \mathbb{R}^n,$$

is locally asymptotically stable if and only if

$$\min_i |\arg(\lambda_i)| > \alpha \frac{\pi}{2}, \quad i = 1, 2, \dots, n. \quad (12)$$

**Theorem 4.** [42] Let E be an equilibrium point of the fractional-order nonlinear system

$$D^{\alpha} = f(x), \quad 0 < \alpha < 2.$$

If the eigenvalues of the Jacobian matrix $A = \frac{\partial f}{\partial x}|_E$ satisfy

$$\min_i |\arg(\lambda_i)| > \alpha \frac{\pi}{2}, \quad i = 1, 2, \dots, n,$$

then the system is asymptotically stable at the equilibrium point E.

It can be verified that the Jacobian matrix of system (10) at the equilibrium point E is

$$J_E = \begin{pmatrix} 1 - x_e^2(b) & -1 \\ \varepsilon & -b\varepsilon \end{pmatrix}$$

and the characteristic polynomial reads

$$P(\lambda) = \lambda^2 + (x_e^2(b) + b\varepsilon - 1)\lambda + \varepsilon(1 + b(x_e^2(b) - 1)).$$

Then using proposition 1, in [47], the fixed point E is locally asymptotically stable if and only if

$$\varepsilon(1 + b(x_e^2(b) - 1)) > 0 \text{ and } (x_e^2(b) + b\varepsilon - 1) > -2\sqrt{\varepsilon(1 + b(x_e^2(b) - 1))\cos(\alpha \frac{\pi}{2})}. \quad (13)$$

Particularly for $(x_e^2(b) - b\varepsilon - 1)^2 < 4\varepsilon$, the Jacobian matrix $J_E$ has a pair of complex conjugate eigenvalues:

$$\lambda_{\pm} = \frac{-(x_e^2(b) + b\varepsilon - 1) \pm i\sqrt{-(x_e^2(b) - b\varepsilon - 1)^2 + 4\varepsilon}}{2}.$$

and the point E is a focus.
## 5 Hopf-Like Bifurcation

An important feature of fractional-order autonomous differential systems is that periodic solution cannot exist [13, 43–45]. Therefore, the classical Hopf bifurcation cannot occur in such systems.

We introduce a novel kind of “bifurcation”, which captures the essential of the Hopf bifurcation; that is the change from static to near periodic solution. We call it Hopf-Like Bifurcation (HLB). The HLB can be characterized when a fixed point of the underlying dynamical system changes its stability property as a pair of complex conjugate eigenvalues $\lambda_{\mp}$ of the Jacobian matrix at the fixed point cross the boundary of an angular sector $|\arg(\lambda_{\mp})| = \alpha\frac{\pi}{2}$ of the complex plane, giving rise (or vanishing) to a small-amplitude S-asymptotically T-periodic solution. Few years ago, some criteria of HLB in fractional-order systems were already introduced by Abdelouahab et al. [46], although it was not called “Hopf-like” therein.

It had been shown recently that Eq. (10) undergoes an HLB at its unique fixed point $E = (x_e(b), y_e(b))$ with respect to the parameter $b$ and the parameter $\alpha$. Of course $E$ is also depending on parameters $a$ and $I$. However, in this study, we fix these parameters and vary only parameters $b$ and $\alpha$.

In order to analyze HLB near this point $E$, we define a function

$$M(b, \alpha) = \alpha \frac{\pi}{2} - \left| \arctan \left( \frac{\sqrt{-(x_e^2(b) - b\varepsilon - 1)^2 + 4\varepsilon}}{x_e^2(b) + b\varepsilon - 1} \right) \right|$$

Then, considering the solution $(b^*, \alpha^*)$ of $M(b, \alpha) = 0$, it is possible to prove the existence of HLB.

**Theorem 5.** (Abdelouahab et al., [14]) Let $(b^*, \alpha^*)$ be a solution to $M(b, \alpha) = 0$. If

$$(x_e^2(b) - b\varepsilon - 1)^2 < 4\varepsilon,$$

and

$$(2x_e(b) \frac{dx_e(b)}{db}(b^2\varepsilon - b(x_e^2(b) - 1) - 2) + (x_e^2(b) - 1)^2 - b\varepsilon(x_e^2(b) - 1) - 2\varepsilon) \Big|_{b=b^*} \neq 0,$$

then system (10) undergoes an HLB at the unique equilibrium point $E$, when $(b, \alpha) = (b^*, \alpha^*)$.

**Numerical example:** Taking the following values of the parameters: $a = 0.75$, $I = 0.41$, $\alpha = 0.05$, $\alpha \in (0, 2)$, $0 < b < 1.4377$, Fig. 3 displays the critical curve $\gamma$ of the following equation:

$$M(b, \alpha) = \alpha \frac{\pi}{2} - \min_i |\arg(\lambda_i(b))| = 0,$$

which separates stable and unstable regions in the $(b, \alpha)$ parameter space. All conditions for HLB [46] are satisfied at each point in the curve $\gamma$, implying that when parameters move from stable to unstable regions in the $(b, \alpha)$ parameter space, the fixed point $E$ loses its stability near the critical curve $\gamma$. This gives rise to small-amplitude oscillatory behavior and allows the possibility of developing fractional-order canard solutions.

## 6 Numerical Analysis of Complex MMO Based on Canard Explosion

In this section, we apply the theory of singularly perturbed systems to investigate the canard phenomenon in the FFHN model and highlight some special kind of bifurcation of MMO versus both parameters $b$ and $\alpha$.

### 6.1 Canard cycles

As $\varepsilon \to 0$ the reduced fractional-order equation of system (10) is

$$\begin{cases} x - \frac{1}{3}x^3 - y + I = 0, \\ D^{\alpha}y = x + a - by. \end{cases} \quad (14)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 3 HLB curve in the (b, α) parameter space.</figcaption></figure>

This equation characterizes the slow dynamics. The fast dynamics is characterized by the fractional-order layer equation

$$ \begin{cases} D^{\alpha}x = x - \frac{1}{3}x^3 + I - y = f(x) - y, \\ D^{\alpha}y = 0. \end{cases} \qquad (15) $$

And the critical manifold is given by

$$ S_0 = \{(x,y) \in \mathbb{R}^2 | y = x - \frac{1}{3}x^3 + I = f(x)\}. \qquad (16) $$

The derivative is $f'(x) = 1 - x^2$, therefore $S_0$ is split into three branches by the two folds $(x,y) = (\pm 1, \pm \frac{2}{3} + I)$, two attractive branches $S_a$, where $f'(x) < 0$, and one repulsive branch $S_r$, where $f'(x) > 0$, as shown in Fig. 4.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Fractional-order fast and slow subsystems of system (10). Single arrows indicate slow motions along the slow curve $S_0$. Double arrows indicate fast motions outside $S_0$, which possesses two attracting branches, $S_a$, and one repelling branch, $S_r$, separated by fold points (red) of the slow curve, corresponding to saddle-node bifurcation points of the fast subsystem.</figcaption></figure>
## 6.2 MMO Based on canard explosion versus both parameters b and α

In this section, we fix $a = 0.75$, $I = 0.41$, $\epsilon = 0.05$. To illustrate the complex canard explosion versus both parameters b and $\alpha$, we study the solution of Eq. (10), when they vary in the rectangle

$$R = \{(b, \alpha) \in \mathbb{R}^2 | 0 \le b \le 1.3 \text{ and } 0.2 \le \alpha \le 1.4\},$$

following the relationship

$$\alpha = -\frac{b}{2} + 1.35, \qquad (17)$$

as illustrated in Fig. 5. HLB point occurs at $A = (b^*, \alpha^*) \approx (0.8183, 0.9409)$, while oscillations can be observed for $b < 0.8183$ and $\alpha > 0.9409$.

The system (10) is numerically integrated on the time interval $[0, t_f]$ using the Grünwald-Letnikov approximation. The step size of this integration is $h = 0.01$ and the initial condition is $(x_0 = x_e(b); y_0 = y_e(b) - 0.005)$. To integrate a fractional derivative system is time consuming; therefore, in order to accelerate the computation the short memory length principle is used, with memory length $L = 100$ [34, 48].

In order to illustrate the complex canard explosion we explore the neighbourhood of the HLB, particularly its-left-hand side segment $[AB]$, where $B \approx (0.7780.961)$ (brown color Fig. 5).

<figure><img src="image_7.png" /><figcaption>FIGURE 5 HLB curve in the $(b, \alpha)$ parameter space.</figcaption></figure>

For $(b, \alpha) = B$, one can only observe large-amplitude oscillations. For $(b, \alpha) = A$, one can only observe small-amplitude oscillations (with amplitude close to zero). When parameters $(b, \alpha)$ are varied from B to A, the number of small-amplitude oscillations, $NSAO((b, \alpha))$, which occurs between every two successive large-amplitude oscillations, changes from $NSAO(B) = 0$ to $NSAO(A) = +\infty$.

To localize infinitesimal sub-segments for which $NSAO((b, \alpha))$ increases, where canard cycles can be developed, the 'Global Local Canard Explosion Search Algorithm' (GLCESA) introduced in [14] is applied. The algorithm is adopted to calculate the number of small-amplitude oscillations belonging between the first and the second large-amplitude oscillations. It determines an infinitesimal parameter sub-segment on which this number increases. A total of 15 canard explosion parameter sub-segments, $CEPS = [\bar{B}_i \bar{A}_i]$, $i = 1, 2, ..., 15$, where $\bar{B}_i = (\bar{b}_i, \bar{\alpha}_i)$ and $\bar{A}_i = (\bar{b}_i + 2 \times 10^{-13}, \bar{\alpha}_i + 10^{-13})$, are displayed in the Table 1, with their corresponding $NSAO$, $tf$ and $PSD$.

From this table and Fig. 6, one can see that the amplitude of the last small oscillation increases, as the parameter $b$ decreases and the parameter $\alpha$ increases. Then, one can observe canard explosion within an exponentially small neighbourhood of each $(\bar{b}_i, \bar{\alpha}_i)$, where the transition from small-amplitude oscillation (stationary-like
<table><caption>TABLE 1 Some canard explosion parameter sub-segments: $CEPS = [(\bar{b}_i, \bar{\alpha}_i) (\bar{b}_i + 2 \times 10^{-13}, \bar{\alpha}_i + 10^{-13})]$, $i = 1, 2, \dots, 13$, with their corresponding NSAO, and $tf$, determined using GLCESA as both parameters $b$ and $\alpha$ are varied.</caption><thead><tr><th>NSAO(α, b)</th><th>tf(α, b)</th><th>α̃<sub>i</sub></th><th>b̃<sub>i</sub></th></tr></thead><tbody><tr><td>14</td><td>702.59</td><td>0.9460520040915</td><td>0.8078959918168</td></tr><tr><td>13</td><td>669.43</td><td>0.9462002574928</td><td>0.8075994850142</td></tr><tr><td>12</td><td>637.49</td><td>0.9463712718829</td><td>0.8072574562340</td></tr><tr><td>11</td><td>609.12</td><td>0.9465701755015</td><td>0.8068596489968</td></tr><tr><td>10</td><td>955.08</td><td>0.9468114549302</td><td>0.8063770901394</td></tr><tr><td>9</td><td>500.67</td><td>0.9470807229186</td><td>0.8058385541626</td></tr><tr><td>8</td><td>469.37</td><td>0.9474164241791</td><td>0.8051671516416</td></tr><tr><td>7</td><td>436.78</td><td>0.9478315727039</td><td>0.8043368545920</td></tr><tr><td>6</td><td>365.56</td><td>0.9483577759722</td><td>0.8032844480554</td></tr><tr><td>5</td><td>332.4</td><td>0.9490323331705</td><td>0.8019353336588</td></tr><tr><td>4</td><td>300.35</td><td>0.9499488046301</td><td>0.8001023907396</td></tr><tr><td>3</td><td>227.11</td><td>0.9512805068416</td><td>0.7974389863166</td></tr><tr><td>2</td><td>189.64</td><td>0.9532285469108</td><td>0.7935429061782</td></tr><tr><td>1</td><td>156.56</td><td>0.9564716090280</td><td>0.7870567819438</td></tr><tr><td>0</td><td>204.82</td><td>0.9608101829226</td><td>0.7783796341546</td></tr></tbody></table>

behavior) to large-amplitude oscillation (relaxation oscillation) happens via a fractional canard solution, as illustrated for example in Fig. 6, where the canard explosion occurs for $(b, \alpha)$ in the parameter sub-segment $[\bar{B}_4, \bar{A}_4]$ with $\bar{B}_4 = (0.7974389863166, 0.9512805068416)$ and $\bar{A}_4 = (0.7974389863168, 0.9512805068415)$.

For $(b, \alpha) = \bar{B}_4$, there are only 3 small-amplitude oscillations between the first and the second large-amplitude oscillations, but for $(b, \alpha) = \bar{A}_4$, there are 5 small-amplitude oscillations.

When a trajectory of a smooth two-dimensional autonomous fractional-order system like FFHN, self-crossing does not imply that the trajectory follows a periodic orbit, because such type of systems cannot have periodic orbit due to the memory effect. This implies that the Poincaré-Bendixson theorem, which prohibits the existence of chaotic orbit in two-dimensional autonomous systems, is only valid for the integer-order setting. When fractional-order dynamical systems are considered, there is more flexibility for the trajectories with the possibility of having chaotic solutions in dimension two, if the fractional-order is close to one. This allows us to the following conjecture:

**Conjecture 1.** *Chaos can exist in two-dimensional autonomous fractional-order dynamical systems like FFHN with the fractional-order close to one.*

## 7 Conclusions

This article has investigated the fractional-order FFHN model to find MMO based on canard explosion when fractional-order is close to one. The study has highlighted the appearance of patterns of solutions with an increasing number of small-amplitude oscillations in each of such patterns, when two parameters vary following a peculiar relationship. To rigorously analyze such complex dynamics of the FFHN model, a new mathematical tool, the Hopf-like bifurcation, is used, which gives a precise definition of the change between a fixed point and an *S*-asymptotically *T*-periodic solution of a fractional-order system. Then, the stability and HLB of the FFHN model have been carefully studied, using a recently introduced algorithm named the Global-Local Canard Explosion Search Algorithm. This analysis confirms the existence of canard oscillations in the neighborhood of a HLB point versus both parameters $b$ and $\alpha$ (the fractional-order). Finally, it was conjectured that chaos can exist in FFHN, which is a very simple model, based on a two-dimensional fractional-order system. It can offer
<figure><img src="image_8.png" /><figcaption>FIGURE 6 Canard solutions observed from the fractional-order system (10): (a) Phase portrait for (b, α) = (0.7974389863166, 0.9512805068416). (b) Time evolution of x for (b, α) = (0.7974389863166, 0.9512805068416). (c) Phase portrait for (b, α) = (0.7974389863168, 0.9512805068415). (d) Time evolution of x for (b, α) = (0.7974389863168, 0.9512805068415).</figcaption></figure>
an incredible possibility for describing the complex dynamics of a neuron, although it is only 2 dimensional.

Références

[1] Hodgkin, A.L., Huxley, A.F. (1952), A quantitative description of membrane current and its application to conduction and excitation in nerve, J. Physiol. 117, 500-544.

[2] Forrest, M. D. (2014), Can the thermodynamic Hodgkin-Huxley model of voltage-dependent conductance extrapolate for temperature ?, Computation, 2, 47-60.

[3] Shilnikov, A. L., Rulkov, N.F. (2003), Origin of chaos in a two-dimensional map modelling spiking-bursting neural activity, International Journal of Bifurcation and Chaos, 13, 3325-3340.

[4] Shilnikov, A. L., Kolomiets, M. (2008), Methods of the qualitative theory for the Hindmarsh-Model: A case study. A Tutorial, International Journal of Bifurcation and Chaos, 18, 2141-2168.

[5] Fitzhugh, R. (1961), Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. 1, 445-466.

[6] Hindmarsh, J.L., Rose, R.M. (1982), A model of the nerve impulse using two first-order differential equations, Nature, 296, 162-164.

[7] Hindmarsh, J.L., Rose, R.M. (1984), A model of neuronal bursting using three coupled first order differential equations, Proc.R.Soc. London B. 221, 87-102.

[8] Perko, L. (2002), Differential Equations and Dynamical Systems, Springer, New York.

[9] Liu, Y., Xie, Y., Kang, Y., Tan, N., Jiang, J., Xu, J.X. (2010), Dynamical characteristics of the fractional-order Fitzhugh-Nagumo model neuron, Advances in Cognitive Neurodynamics (II), In: Wang R., Gu Fanji (eds.), pp. 253-258, Springer, Dordrecht (2010).

[10] Momani, S. A. Freihat, M. Al-Smadi (2014), Analytical study of fractional-order multiple chaotic Fitzhugh-Nagumo neurons model using multistep generalized differential transform method, Abstract Appl. Anal. 24, ID 276279.

[11] Armanyos, M., Radwan, A.G. (2016), Fractional-order Fitzhugh-Nagumo and Izhikevich neuron models, Proc. 13th Int. Conf. on Elect. Engin./Electron., Computer, Telecomm. Info. Tech. (ECTI-CON), pp. 1-5.

[12] Alidousti, J., Khoshsiar Ghaziani, R. (2017), Spiking and bursting of a fractional-order of the modified Fitzhugh-Nagumo neuron model, Math. Models Comp. Simul. 9(3), 390-403.

[13] Tavazoei, M. S., Haeri, M. (2009), A proof for non-existence of periodic solutions in time invariant fractional-order systems, Automatica, 45, 1886-1890.

[14] Abdelouahab M-S, Lozi R, Chen G. (2019), Complex canard explosion in a fractional-order Fitzhugh-Nagumo model. Int. J. Bifurc. Chaos, 29(8): 1950111-1950133.

[15] Zhabotinsky, A. M., (1964), Periodic kinetics of oxidation of malonic acid in solution, Biofizika, 9, 306-311.

[16] Brøns, M., Krupa, M., Wechselberger, M., (2006), Mixed mode oscillations due to the ge,realized canard phenomenon, Fields Institute Communications, 49, 39-63.

[17] Cartier, P., (1982), Perturbations singulières des équations différentielles ordinaires et analyse non-standard, Astérisque, tome 92-93, Séminaire Bourbaki, exp. no 580, p. 21-44.

[18] Benoît, E., Callot, J. F., Diener, F., Diener, M. (1981), Chasse au canard, Collect. Math. 31, 37-119.

[19] Desroches, M., Jeffrey, M. R. (2011), Canards and curvature: The smallness of $\epsilon$ in slow-fast dynamics, Proc. Roy. Soc., A 467, 2404-2421.

[20] Van der Pol, B. (1926), On relaxation oscillations, London, Edinburgh, and Dublin Phil. Mag. J. Sci. 7, 978-992.

[21] Wechselberger, M. (2012), À propos de canards (A propos canards), Trans. Am. Math. Soc. 364, 3289-3309.

[22] Shchepakina, E., Sobolev, V., Mortell, M.P. (2014), Singular Perturbations: Introduction to System Order Reduction Methods with Applications, Lecture Notes in Mathematics, Springer.

[23] Ermentrout, G. B., Wechselberger, M. (2009), Canards, clusters, and synchronization in a weakly coupled interneuron model, SIAM J. Appl. Dyn. Syst., 8, 253-278.

[24] Rubin, J., Wechselberger, M. (2008), The selection of mixed-mode oscillations in a Hodgkin-Huxley model with multiple timescales, Chaos, 18, 015105.

[25] Guckenheimer, J., Scheper, C. (2013), Multiple time scale analysis of a model Belousov-Zhabotinskii reaction, SIAM J. Appl. Dyn. Syst., 12, 1968-1996.

[26] Marino, F., Ciszak, M., Abdalah, S. F., Al-Naimee, K., Meucci, R., Arecchi, F. T. (2011), Mixed-mode oscillations via canard explosions in light-emitting diodes with optoelectronic feedback, Phys. Rev. E 84, 047201.

[27] Rinzel, J., Miller, R. (1980), Numerical calculation of stable and unstable periodic solutions to the Hodgkin-Huxley equations, Math. Biosci., 49, 27-59.

[28] Rubin, J., Wechselberger, M. (2007), Giant squid-hidden canard: The 3D geometry of the Hodgkin Huxley model, Biol. Cybern., 97, 5-32.

[29] Leibniz, G.W. (1626), Leibnizens mathematische Schriften, Georg Ohms Verlagsbuch Handlung Hildesheim.

[30] Bagley, R. L., Calico, R. A. (1991), Fractional-Order state equations for the control of viscoelastically damped struc-
tures, J. Guid. Contr. Dyn., 14, 304-311.

[31] Ichise, M., Nagayanagi, Y., Kojima, T. (1971), An analog simulation of noninteger order transfer functions for analysis of electrode process, J. Electro. Chem., 33, 253-65.

[32] Heaviside, O. (1971), Electromagnetic Theory, Chelsea, New York.

[33] Kusnezov, D., Bulgac, A., Dang, G. D. (1999), Quantum levy processes and fractional kinetics, Phys. Rev. Lett., 82, 1136-1139.

[34] Podlubny, I. (1999), Fractional Differential Equations Academic Press, San Diego.

[35] Samko, S.G., Kilbas, A.A., Marichev, O.I. (1993), Fractional integrals and derivatives: Theory and Applications, Gordan and Breach, Amsterdam.

[36] Butzer, P.L., Westphal, U. (2000), An introduction to fractional calculus, In: R. Hilfer (ed.), Applications of Fractional Calculus in Physics, World Scientific, Singapore, pp. 1-85.

[37] Caputo, M. (1967), Linear models of dissipation whose Q is almost frequency independent-II, Geophys J. R. Astron. Soc. 13, 529-539.

[38] Nagumo, J., Arimoto, S., Yoshizawa, S. (1962), An active pulse transmission line simulating nerve axon, Proc. IRE 50, 2061-2070.

[39] Jonscher, A. K. (1983), Dielectric Relaxation in Solids Chelsea Dielectric Press, London.

[40] Westerlund, S., Ekstam, S. (1994), Capacitor theory, IEEE Trans. Dielectr. Electr. Insul. 1, 826-839.

[41] Westerlund, S. (1991), Dead matter has memory!, Physica Scripta, 43, 174.

[42] Abdelouahab, M-S., Hamri, N.E., Wang, J. (2010), Chaos control of a fractional-order financial system, Math. Probl. in Eng. ID 270646, 1-18.

[43] Tavazoei, M-S. (2010), A note on fractional-order derivatives of periodic functions, Automatica, 46, 945-948.

[44] Kaslik, E., Sivasundaram, S. (2012), Non-existence of periodic solutions in fractional-order dynamical systems and a remarkable difference between integer and fractional-order derivatives of periodic functions, Nonl. Anal.: Real World Appl. 13, 1489-1497.

[45] Kang, Y.M., Xie, Y., Lu, J. C., Jiang, J. (2015), On the non-existence of non-constant exact periodic solutions in a class of the Caputo fractional-order dynamical systems, Nonlin. Dyn. 82, 1259-1267.

[46] Abdelouahab, M-S., Hamri, N. E., Wang, J. W. (2012), Hopf bifurcation and chaos in fractional-order modified hybrid optical system, Nonlin. Dyn. 69, 275-284.

[47] Bourafa, S., Abdelouahab, M-S., Moussaoui A. (2020), On some extended Routh-Hurwitz conditions for fractional-order autonomous systems of order $\alpha \in [0,2)$ and their applications to some population dynamic models, Chaos, Solitons and Fractals, 109623.

[48] Abdelouahab, M. S., Hamri, N. E. (2016), The Grünwald-Letnikov fractional-order derivative with fixed memory length, Mediterr. J. Math. 13, 557-572.

[49] Matignon, D. (1996), Stability results in fractional differential equation with applications to control processing, Proc. multiconf. comput. Eng. in Syst. Appl. (IMICS), IEEE-SMC. 2, 963-968.

[50] Moze, M., Sabatier, J. (2005), LMI tools for stability analysis of fractional systems, Proc. ASME Int. Design Eng. Tech. Conf. & Comput. Infor. in Eng. Conf., Long Beach, CA, 24-28 Sept.