This work is licensed under the Creative Commons Attribution 4.0 International License.
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, 10, 11, 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. 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.
One 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.
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.
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
\left\{{\matrix{{\dot x = y - {1 \over 3}{x^3} + x{\kern 1pt},} \cr {\dot y = \varepsilon (a - x){\kern 1pt}.\;\;\;\;\;} \cr}} \right.
For a small positive parameter ɛ ≪ 1, the variable x is driven by the fast vector field evolving on the fast time scale t = τ/ɛ, and y evolves on a slow time scale τ directed by the slow vector field. Thus, x (resp. y) is called the fast (resp. slow) variable. The cubic shaped curve S, defined by ẋ = 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].
Figure 1
Canard explosion of the van der Pol oscillator for ɛ = 0.02 happens in an exponentially small parameter interval near a = a* ≈ 0.997461066999 where the transition from relaxation oscillations for (a) a < 0.997461066999, to small-amplitude limit cycles for (c) a ≥ 0.997461066999, comes via (b) canard cycles.
In general, canard solutions occur in singularly perturbed systems of the form
\left\{{\matrix{{\varepsilon {{dx(\tau)} \over {d\tau}} = f(x(\tau),y(\tau),\varepsilon){\kern 1pt},} \cr {{{dy(\tau)} \over {d\tau}} = g(x(\tau),y(\tau),\varepsilon){\kern 1pt},} \cr}} \right.
with x ∈ ℝn and y ∈ ℝm being the fast (resp. the slow) variable, both functions f and g are sufficiently smooth, and τ denotes the independent variable of (2), called the slow time scale. After a time, rescaling to the fast time scale t = τ/ɛ (for ɛ ≠ 0), (2) is equivalent to the system
\left\{{\matrix{{\dot x = f(x,y,\varepsilon){\kern 1pt},} \cr {\dot y = \varepsilon g(x,y,\varepsilon){\kern 1pt},} \cr}} \right.
where the dot denotes differentiation with respect to t.
As, ɛ → 0, system (2) is reduced to the slow subsystem
\left\{{\matrix{{0 = f(x,y,0){\kern 1pt},} \cr {\dot y = g(x,y,0){\kern 1pt},} \cr}} \right.
and system (3) to the fast (or layer) subsystem
\left\{{\matrix{{\dot x = f(x,y,0),} \hfill \cr {\dot y = 0.} \hfill \cr}} \right.
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 ɛ > 0 [21]. The critical (slow) manifold
S = \left\{{(x,y) \in \mathbb {R}^n \times \mathbb {R}^m:f(x,y,0) = 0} \right\}
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 Sa ⊂ S, for which all eigenvalues of Dxf have negative real parts, is called the attracting slow invariant manifold and, similarly, the subset Sr ⊂ S for which all eigenvalues of Dxf 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 Sa and then continuing for a while along Sr, is called a canard solution. Conversely, the solution first follows Sr, and then continuing after Sr for a while, is called a false canard solution.
Mixed-mode oscillations
MMO consists of L large-amplitude (relaxation) oscillations followed by s small-amplitude (subthreshold) oscillations, simply denoted by Ls [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].
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, 31, 32, 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 Fitzugh-Nagumo system [5], with only two equations.
Several definitions of fractional-order derivatives have been introduced in the literature [34, 35, 36, 37]. A common one is the Riemann-Liouville definition of fractional derivatives [34, 35], given by
\matrix{{\ _a^RD_t^\alpha x(t)} \hfill & {= {1 \over {\Gamma (m - \alpha)}}{{{d^m}} \over {d{t^m}}}\int_a^t {{(t - \tau)}^{m - \alpha - 1}}x(\tau)d\tau} \hfill \cr {} \hfill & {= {{{d^m}} \over {d{t^m}}}{(_a}j_t^{m - \alpha}x(t)),\;\;t > a,\;m - 1 \le \alpha < m{\kern 1pt},} \hfill \cr}
where Γ is the gamma function and
{\ _a}j_t^\beta
is the Riemann-Liouville integral operator defined as
{\ _a}j_t^\beta x(t) = {1 \over {\Gamma (\beta)}}\int_a^t {(t - \tau)^{\beta - 1}}x(\tau)d\tau {\kern 1pt}.
The Laplace transform of the α-order Riemann–Liouville differential operator is
L\left\{{\ _0^RD_t^\alpha x(t)} \right\} = {s^\alpha}L\left\{{x(t)} \right\} - \sum\limits_{k = 0}^{m - 1} {s^k}{\left[ {\ _0^RD_t^{\alpha - 1 - k}x(t)} \right]_{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:
\matrix{{\ _a^cD_t^\alpha x(t)} \hfill & {= {1 \over {\Gamma (n - \alpha)}}\int_a^t {{(t - \tau)}^{n - \alpha - 1}}{x^{(n)}}(\tau)d\tau} \hfill \cr {} \hfill & {= {j^{n - \alpha}}\left({{{{d^n}} \over {d{t^n}}}x(t)} \right),\;t > a{\kern 1pt},} \hfill \cr}
where n = ⌈α⌉ is the value of α rounded up to the nearest integer. The Laplace transform of the α-order Caputo differential operator is
L\left\{{\ _0^CD_t^\alpha x(t)} \right\} = {s^\alpha}L\left\{{x(t)} \right\} - \sum\limits_{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
\ _a^GD_t^\alpha x(t) = \mathop {\lim}\limits_{h \to 0} {1 \over {{h^\alpha}}}\sum\limits_{k = 0}^{k = {{t - a} \over h}} {(- 1)^k}{{\Gamma (\alpha + 1)} \over {k!\Gamma (\alpha - k + 1)}}x(t - kh){\kern 1pt},
where t > a and α is a positive real number. Its integral form is
\ _a^GD_t^{- \alpha}x(t) = \mathop {\lim}\limits_{h \to 0} {h^\alpha}\sum\limits_{k = 0}^{k = {{t - a} \over h}} {{\Gamma (\alpha + k)} \over {k!\Gamma (\alpha)}}x(t - kh){\kern 1pt}.
When x is of class Cm, where m − 1 ≤ α < 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.
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.
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.
Figure 2
Electrical circuit equivalent to the Fitzhugh–Nagumo model.
The equation of this model reads
\left\{{\matrix{{{{dv} \over {dt}} = v - {{{v^3}} \over 3} - w + I} \cr {{{dw} \over {dt}} = {1 \over T}(v + a - bw)} \cr}} \right.
with parameters I, a, b and T.
On can reduce it to the form of Eq. (1), by posing x = v, y = w and ɛ = 1/T, as
\left\{{\matrix{{\dot x = x - {{{x^3}} \over 3} - y + I,} \cr {\dot y = \varepsilon (x + a - by).} \cr}} \right.
It is obvious that, if b = I = 0 and x → −x, Eq. (9) is exactly the van der Pol system (1).
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{{dv(t)} \over {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) = C{D^\alpha}v(t),
where α 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) = L{D^\alpha}I(t),
where the derivative order α 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):
\left\{{\matrix{{{D^{{\alpha _1}}}x = x - {1 \over 3}{x^3} - y + I,} \hfill \cr {{D^{{\alpha _2}}}y = \varepsilon (x + a - by),} \hfill \cr}} \right.
where the derivative order α1 (resp. α2 is constant and related to the loss of the capacitor (resp. the proximity effect of the inductor).
Stability of the fixed point of FFHN
In this section, we briefly recall some results concerning the fixed points of the FFHN for α1 = α2 = α ∈ (0, 2). In order to have a unique fixed point, some parameter restrictions are introduced.
[14] The fixed point E of Eq. (10) is locally asymptotically stable if\left| {\arctan \left({{{\sqrt {- {{(x_e^2(b) - b\varepsilon - 1)}^2} + 4\varepsilon}} \over {x_e^2(b) + b\varepsilon - 1}}} \right)} \right| > \alpha {\pi \over 2}{\kern 1pt}.
The proof is straightforward based on the following theorems.
Theorem 3
[49, 50] The following fractional-order linear autonomous system:
\left\{{\matrix{{{D^\alpha}X = AX{\kern 1pt},} \hfill \cr {X(0) = {X_0}{\kern 1pt},} \hfill \cr}} \right.\;\;X \in {{\bf{R}}^n},\;0 < \alpha < 2\;\;and\;A \in {{\bf{R}}^n} \times {{\bf{R}}^n},is locally asymptotically stable if and only if\mathop {\min}\limits_i \left| {\arg ({\lambda _i})} \right| > \alpha {\pi \over 2},\;i = 1,2,...,n{\kern 1pt}.
Theorem 4
[42] Let E be an equilibrium point of the fractional-order nonlinear system{D^\alpha} = f(x){\kern 1pt},\;\;0 < \alpha < 2{\kern 1pt}.
If the eigenvalues of the Jacobian matrixA = {\left. {{{\partial f} \over {\partial x}}} \right|_E}satisfy\mathop {\min}\limits_i \left| {\arg ({\lambda _i})} \right| > \alpha {\pi \over 2},\;i = 1,2,...,n{\kern 1pt},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} = \left({\matrix{{1 - x_e^2(b)} & {- 1} \cr \varepsilon & {- b\varepsilon} \cr}} \right){\kern 1pt}
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)){\kern 1pt}.
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\;{\rm{and}}\;(x_e^2(b) + b\varepsilon - 1) > - 2\sqrt {\varepsilon (1 + b(x_e^2(b) - 1))} \cos (\alpha {\pi \over 2}).
Particularly for
{(x_e^2(b) - b\varepsilon - 1)^2} < 4\varepsilon
, the Jacobian matrix JE has a pair of complex conjugate eigenvalues:
{\lambda _ \pm} = {{- (x_e^2(b) + b\varepsilon - 1) \pm i\sqrt {- {{(x_e^2(b) - b\varepsilon - 1)}^2} + 4\varepsilon}} \over 2}{\kern 1pt}.
and the point E is a focus.
Hopf-Like Bifurcation
An important feature of fractional-order autonomous differential systems is that periodic solution cannot exist [13, 43, 44, 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 λ∓ of the Jacobian matrix at the fixed point cross the boundary of an angular sector
\left| {\arg ({\lambda _{\; \mp}})} \right| = \alpha {\pi \over 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 = (xe(b), ye(b)) with respect to the parameter b and the parameter α. E is also depending on parameters a and I. However, in this study, we fix these parameters and vary only parameters b and α.
In order to analyze HLB near this point E, we define a function
M(b,\alpha) = \alpha {\pi \over 2} - \left| {\arctan \left({{{\sqrt {- {{(x_e^2(b) - b\varepsilon - 1)}^2} + 4\varepsilon}} \over {x_e^2(b) + b\varepsilon - 1}}} \right)} \right|{\kern 1pt}.
Then, considering the solution (b*, α*) of M(b, α) = 0, it is possible to prove the existence of HLB.
Theorem 5
(Abdelouahab et al., [14]) Let (b*, α*) be a solution to M(b, α) = 0. If{(x_e^2(b) - b\varepsilon - 1)^2} < 4\varepsilon,and{\left. {\left({2{x_e}(b){{d{x_e}(b)} \over {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} \right)} \right|_{b = {b^ *}}} \ne 0,then system (10) undergoes an HLB at the unique equilibrium point E, when (b, α) = (b*, α*).
Numerical example: Taking the following values of the parameters: a = 0.75, I = 0.41, α = 0.05, α ∈ (0, 2), 0 < b < 1.4377, Fig. 3 displays the critical curve γ of the following equation:
M(b,\alpha) = \alpha {\pi \over 2} - \mathop {\min}\limits_i \left| {\arg ({\lambda _i}(b))} \right| = 0{\kern 1pt},
which separates stable and unstable regions in the (b, α) parameter space. All conditions for HLB [46] are satisfied at each point in the curve γ, implying that when parameters move from stable to unstable regions in the (b, α) parameter space, the fixed point E loses its stability near the critical curve γ. This gives rise to small-amplitude oscillatory behavior and allows the possibility of developing fractional-order canard solutions.
Figure 3
HLB curve in the (b, α) parameter space.
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 α.
Canard cycles
As ɛ → 0 the reduced fractional-order equation of system (10) is
\left\{{\matrix{{x - {1 \over 3}{x^3} - y + I = 0{\kern 1pt},} \cr {{D^\alpha}y = x + a - by.\;\;} \cr}} \right.
This equation characterizes the slow dynamics. The fast dynamics is characterized by the fractional-order layer equation
\left\{{\matrix{{{D^\alpha}x = x - {1 \over 3}{x^3} + I - y = f(x) - y{\kern 1pt},} \hfill \cr {{D^\alpha}y = 0{\kern 1pt}.} \hfill \cr}} \right.
And the critical manifold is given by
{S_0} = \{(x,y) \in \mathbb {R}^2|\;y = x - {1 \over 3}{x^3} + I = f(x)\} {\kern 1pt}.
The derivative is f′(x) = 1 − x2, therefore S0 is split into three branches by the two folds
(x,y) = (\pm 1, \pm {2 \over 3} + I)
, two attractive branches Sa, where f′(x) < 0, and one repulsive branch sr, where f′(x) > 0, as shown in Fig. 4.
Figure 4
Fractional-order fast and slow subsystems of system (10). Single arrows indicate slow motions along the slow curve S0. Double arrows indicate fast motions outside S0, which possesses two attracting branches, Sa, and one repelling branch, Sr, separated by fold points (red) of the slow curve, corresponding to saddle-node bifurcation points of the fast subsystem.
MMO Based on canard explosion versus both parameters b and α
In this section, we fix a = 0.75, I = 0.41, ɛ = 0.05. To illustrate the complex canard explosion versus both parameters b and α, we study the solution of Eq. (10), when they vary in the rectangle
R = \left\{{(b,\alpha) \in \mathbb {R}^2|\;0 \le b \le 1.3\;and\;0.2 \le \alpha \le 1.4} \right\},
following the relationship
\alpha = - {b \over 2} + 1.35,
as illustrated in Fig. 5. HLB point occurs at A = (b*, α*)≈ (0.8183, 0.9409), while oscillations can be observed for b < 0.8183 and α > 0.9409.
Figure 5
HLB curve in the (b, α) parameter space.
The system (10) is numerically integrated on the time interval [0, tf] using the Grünwald-Letnikov approximation. The step size of this integration is h = 0.01 and the initial condition is (x0 = xe(b);y0 = ye(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 ≈ (0.7780.961) (brown color Fig. 5).
For (b, α) = B, one can only observe large-amplitude oscillations. For (b, α) = A, one can only observe small-amplitude oscillations (with amplitude close to zero). When parameters (b, α) are varied from B to A, the number of small-amplitude oscillations, NSAO((b, α)), which occurs between every two successive large-amplitude oscillations, changes from NSAO(B) = 0 to NSAO(A) = +∞.
To localize infinitesimal sub-segments for which NSAO((b, α)) 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.
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,..., 13, with their corresponding NSAO, and tf, determined using GLCESA as both parameters b and α are varied.
NSAO(α, b)
t f (α, b)
{\bar \alpha _i}
{\bar b_i}
14
702.59
0.9460520040915
0.8078959918168
13
669.43
0.9462002574928
0.8075994850142
12
637.49
0.9463712718829
0.8072574562340
11
609.12
0.9465701755015
0.8068596489968
10
955.08
0.9468114549302
0.8063770901394
9
500.67
0.9470807229186
0.8058385541626
8
469.37
0.9474164241791
0.8051671516416
7
436.78
0.9478315727039
0.8043368545920
6
365.56
0.9483577759722
0.8032844480554
5
332.4
0.9490323331705
0.8019353336588
4
300.35
0.9499488046301
0.8001023907396
3
227.11
0.9512805068416
0.7974389863166
2
189.64
0.9532285469108
0.7935429061782
1
156.56
0.9564716090280
0.7870567819438
0
204.82
0.9608101829226
0.7783796341546
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 α 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 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, α) 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).
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).
For (b, α) =
{\bar B_4}
, there are only 3 small-amplitude oscillations between the first and the second large-amplitude oscillations, but for (b, α) =
{\bar A_4}
, there are 5 small-amplitude oscillations.
More complex oscillations
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 this type of systems cannot have periodic orbits 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 valid only 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 brings up 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.
Some examples of chaotic solutions are shown for highlighting this conjecture.
More complex oscillations
In the neighborhoods of HLB points, more complex oscillations arise. Successive MMO Ls are organized in a mixture of repetitive patterns LsLs+1 or LsLsLs+1 like 15 − 15 − 16, 15 − 15 − 16, 15 − 16, 15 − 16, 15 − 16, 15 − 16, 15 − 16, 15 − 16 as shown in Fig. 7 for α = 0.9490476218825, b = 0.8019047562350.
Figure 7
Repetitive patterns of MMO 15 − 15 − 16 and 15 − 16 for α = 0.9490476218825, b = 0.8019047562350
Sometimes the mixing between patterns LsLs+1 and LsLs is not so regular like 14 − 13, 14 − 13, 14 − 14, 14 − 13, 14 − 13, 14 − 13, 14 − 13, 14 − 13 for α = 0.9512573981986, b = 0.7974852036028 as displayed in Fig. 8.
Figure 8
Repetitive patterns of MMO 14 − 13 and 14 − 14 for α = 0.9512573981986, b = 0.7974852036028
Chaotic MMO
Depending upon the initial values x0, y0 of Eq. (14), the MMO are different. Even with an incredibly small difference between such initial conditions, one can find, non super-posable MMO, although forming the same pattern (Fig. 9).
Figure 9
Nonidentical MMO 41, 31, 21, 11 for α = 0.9608101829227, b = 0.7783796341546, red (x0, y0) = (−0.94, −0.27), blue (x0, y0) = (−0.94, −0.26).
According to Strogatz [51], chaos requires three ingredients:
a deterministic system,
aperiodic long-term behavior,
sensitive dependence on initial conditions.
The first condition is obviously verified by Eq. (14), while the other conditions seem present from preliminary numerical experiments. In order to confirm the existence of chaotic MMO organized in patterns, it is necessary to compute Lyapunov exponents. There exists very efficient method to perform such computation for mappings [52], however in the case of fractional-order dynamical systems, this kind of method must be developed.
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 α (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 an incredible possibility for describing the complex dynamics of a neuron, although it is only two dimensional.
