Fractional Calculus-Based Generalization of the FitzHugh-Nagumo Model: Biophysical Justification, Dynamical Analysis and Neurocomputational Implications

Abstract

In this chapter, the dynamical behavior of the incommensurate fractional-order FitzHugh-Nagumo model of neuron is explored in details from local stability analysis. First of all, considering that the FitzHugh-Nagumo model is a mathematical simplification of the Hodgkin-Huxley model, the considered model is derived from the fractional-order Hodgkin-Huxley model obtained taking advantage of the powerfulness of fractional derivatives in modeling certain biophysical phenomena as the dielectrics losses in cell membranes, and the anomalous diffusion of particles in ion channels. Then, it is shown that the fractional-order FitzHugh-Nagumo model can be simulated by a simple electrical circuit where the capacitor and the inductor are replaced by corresponding fractional-order electrical elements. Then, the local stability of the model is studied using the Theorem on the stability of incommensurate fractional-order systems combined with the Cauchy’s argument Principle. At last, the dynamical behavior of the model are investigated, which confirms the results of local stability analysis. It is found that the simple model can exhibit, among others, complex mixed mode oscillations, phasic spiking, first spike latency, and spike timing adaptation. As the dynamical richness of a neuron expands its computational capacity, it is thus obvious that the fractional-order FitzHugh-Nagumo model is more computationally efficient than its integer-order counterpart.

1. Introduction

When excited sufficiently, the neuron, which is the primary unit of brain and nerves electrical activity, generates an action potential, also known as neuronal spike, which is a rapid increase then decrease in the neuronal voltage (see for example [1]). This action potential is at the basis of many mechanisms of communication between neurons and therefore is fundamental to understanding signal processing in the brain and nerves activity [2, 3]. Indeed, action potentials can propagate in essentially constant shape away from the cell body along nerves axons and toward synaptic connections with other cells [3]. Neurons exhibit many different spike-based behaviors including regular periodic and chaotic spiking (train of spikes) and bursting (alternation between a quiescent state and repetitive spiking) [3]. Bursting activity via action potentials plays a crucial role in neuronal communication, including robust transmission of signals [4, 5]. Another interesting spike-based behavior is the mixed mode oscillations (MMOs) pattern, which is an alternation between oscillations of distinct large and small amplitudes [6], where large amplitude oscillations are spikes. Certain MMOs patterns are considered as a type of bursting patterns where the period of small amplitude oscillations is considered as the quiescent phase of the bursting pattern.

Since the pioneering work of Hodgkin and Huxley [1], many relevant conductance-based models and simplified phenomenological models have been developed to describe the brain and nerves microscopic dynamical functions. Studying these computational models with the help of the tools of nonlinear dynamics and singular perturbation theory has shed light on the dynamical processes that support the generation of spiking, bursting and MMOs. Indeed, all these models are nonlinear and almost all of them are singularly perturbed systems which are systems involving multiple time scales. The simplest singularly perturbed systems evolve on on two time scales and are therefore named slow-fast systems. Spike as a form of relaxation, bursting and MMOs are the dynamical signatures of the slow-fast property of a system. The complex slow-fast dynamical behaviors that are bursting and MMOs have been intensively investigated in three-dimensional slow-fast systems and have been found previously to occur only in slow-fast systems involving at least three variables because the successive trigger and extinction of spikes suppose the multiple time scale trip of the system trajectory on a complex well organized high-dimensional phase space. A long time afterwards, it was found that noise can induce MMOs in simple two-variable systems as the van der Pol oscillator with constant forcing [7] and FitzHugh-Nagumo-like oscillators [8, 9]. More recently, we found that MMOs can emerge in two-variable systems due to fractional derivation, while studying the fractional-order van der Pol oscillator with constant forcing where the forcing value was set near the Hopf bifurcation point [10]. Subsequently, this result was confirmed by Abdelouahab et al. while studying Hopf-like bifurcation and bifurcating oscillatory states in a fractional-order FitzHugh-Nagumo model [11].

Recent studies have shown that the fractional derivation increases the dynamical richness of neuronal models. For example, in Ref. [12], Teka et al. have implemented the fractional dynamics on each gating variable of the Hodgkin-Huxley model and they found that the obtained fractional-order versions of the Hodgkin-Huxley model exhibit not only spiking behavior as their integer-order counterpart, but in addition, square wave bursting, MMOs, pseudo-plateau potentials, and phasic spiking. In Ref. [13], Shi and Wang considered a fractional-order Morris-Lecar model and found that this new model exhibits a rich variety of bursting patterns that appear in some common neuron models with properly chosen parameters but do not exist in the corresponding integer-order Morris-Lecar model. In Ref. [14], Mondal et al. considered a fractional-order FitzHugh-Rinzel model whose integer-order counterpart is an elliptic burster and found that it exhibits a wide range of neuronal responses including regular spiking, fast spiking, bursting, and MMOs. In Ref. [15], Teka et al. studied the fractional-order Izhikevich model and found that the model produces a wide range of neuronal spike responses, including regular spiking, fast spiking, intrinsic bursting, MMOs, regular bursting and chattering, by adjusting only the fractional derivatives order. The dynamical richness of these fractional-order systems with at least three variables would justify the occurrence of MMOs in the aforementioned two-variable fractional-order slow-fast systems.

In this chapter, we explore further the dynamical behaviors of fractional-order two-variable slow-fast systems by considering an incommensurate fractional-order FitzHugh-Nagumo (FHN) model derived on the basis of biophysical concepts. Commensurate fractional-order FHN models that are widely considered in previous works are just mathematical generalizations that are not supported by any biophysical justification (see for example Ref. [11]). Exploring the behavior of the incommensurate fractional-order FHN model, we have found that depending on the orders of fractional derivatives the model can exhibit two types of MMOs. In the first case, they are identical to classical folded nodes-induced MMOs observed in integer-order systems, also known as canard generated MMOs [16]. In the second case, obtained for lower derivatives orders, the small oscillations of the MMOs pattern start with very low amplitude which then grows slowly before the oscillations enter the spiking phase. This last type of MMOs is identical to singular Hopf bifurcation-induced MMOs observed in integer-order systems [16]. In the second case, the active MMOs phase is sometimes, depending on initial conditions, preceded by a quiescent state which is known in the context of neuroscience as first spike latency. And in addition, the MMOs regime shows a prominent spike timing adaptation as the order of fractional derivatives decrease. The rest of the paper is organized as follows. In Section 2, the incommensurate fractional-order FHN model is derived from the fractional-order Hodgkin-Huxley model obtained using biophysical concepts. In Section 3, it is shown that the fractional-order FHN model can be simulated using an electrical circuit that is similar to the circuit proposed by Nagumo et al. [17] with the only difference that the capacitor and the inductor are replaced by corresponding fractional-order electrical elements. Section 4 is devoted to the study of local stability of the fractional-order FHN model, which allows to derive the conditions of occurrence of Hopf-like bifurcation with respect to fractional derivatives orders. In Section 5, the dynamical behavior of the fractional-order FHN model is revealed in order to confirm the results of local stability analysis carried out in the preceding section. A particular attention is granted to MMOs and first spike latency as they are new dynamical features due to fractional derivatives. The chapter ends with a conclusion where the results are summarized.

2. Model and biophysical justification

The FHN model

$$\begin{array}{c}\epsilon \frac{\mathit{dx}}{\mathit{dt}}=y+x-\frac{1}{3}{x}^3+I,\\ \frac{\mathit{dy}}{\mathit{dt}}=-x-\mathit{\delta y}+\gamma ,\end{array}$$E1

is a simple representative of a class of excitable-oscillatory systems including the Hodgkin-Huxley model of the squid giant axon. The derivation of the FHN model as a simple nerve membrane model is based on phase space methods applied to the Hodgking-Huxley model. Indeed, the phase diagram of the FHN model and a properly chosen projection from the 4-dimensional Hodgking-Huxley phase space onto a plane are similar, which underlines the relationship between the two models.

As the FHN model is just a mathematical representation, to derive its generalization with fractional derivatives, we consider as starting point the generalization of the Hodgking-Huxley model which was derived using biophysical concepts. The Hodgking-Huxley model is given by the following set of differential equations:

$$\begin{array}{c}{C}_m\frac{\mathit{dv}}{\mathit{dt}}=I-{\overline{g}}_{\mathrm{K}}{n}^4\left(v-v_{\mathrm{K}}\right)-{\overline{g}}_{\mathrm{Na}}{m}^{3}h\left(v-v_{\mathrm{Na}}\right)-{\overline{g}}_l\left(v-v_l\right),\\ \frac{\mathit{dn}}{\mathit{dt}}=a_n\left(1-n\right)-b_{n}n,\\ \frac{\mathit{dm}}{\mathit{dt}}=a_m\left(1-m\right)-b_{m}m,\\ \frac{\mathit{dh}}{\mathit{dt}}=a_h\left(1-h\right)-b_{h}h,\end{array}$$E2

where $v=E_m-E_r$ is the displacement of the membrane potential $E_m$ from its resting value $E_r$; $n$, $m$, and $h$ are the potassium current activation, sodium current activation, and sodium current inactivation gating variables, respectively. $C_m$ is the membrane capacity per unit area, $I$ is an applied stimulus current. $v_{\mathrm{K}}=E_{\mathrm{K}}-E_r$, $v_{\mathrm{Na}}=E_{\mathrm{Na}}-E_r$, and $v_l=E_l-E_r$, where $E_{\mathrm{K}}$ and $E_{\mathrm{Na}}$ are the equilibrium potentials for the sodium and potassium ions, and $E_l$ is the potential at which the leakage current due to chloride and other ions is zero. ${\overline{g}}_{\mathrm{K}}{n}^4=g_{\mathrm{K}}$, ${\overline{g}}_{\mathrm{Na}}{m}^{3}h=g_{\mathrm{Na}}$, where $g_{\mathrm{K}}$, $g_{\mathrm{Na}}$ and ${\overline{g}}_l$ are ionic conductances. $a_n$, $b_n$, $a_m$, $b_m$, $a_h$, $b_h$ are rate constants which vary with membrane voltage. However, Hodgking and Huxley claimed that “The only major reservation which must be made about [the first equation in Eq. (2): $I=C_m\frac{\mathit{dv}}{\mathit{dt}}+I_{\mathrm{Na}}+I_{\mathrm{K}}+I_l$] is that it takes no account of dielectric loss in the membrane” [1]. Now, to account for dielectric losses in a capacitor, Curie proposed the following empiric relation between a DC voltage $U_0$ applied at $t=0$ and the current $i$ it will produce [18]:

$$i\left(t\right)=\frac{U_0}{h_1{t}^{\alpha }},$$E3

where $h_1$ is a constant related to the capacitance of the capacitor and the kind of dielectric, $0<\alpha <1$, $\alpha$ is related to the losses of the capacitor. The lower the losses, the closer to 1 is $\alpha$. Westerlund and Ekstam [19] showed that for a general input voltage $v\left(t\right)$, Eq. (3) turns into:

$$i\left(t\right)=C_{\alpha }\frac{d^{\alpha }v\left(t\right)}{{\mathit{dt}}^{\alpha }},$$E4

where $C_{\alpha }=\mathrm{\Gamma }\left(1-\alpha \right)/h_1$, $\frac{d^{\alpha }v\left(t\right)}{{\mathit{dt}}^{\alpha }}$ is the fractional ($\alpha$-order) time derivative of $v\left(t\right)$, and $\mathrm{\Gamma }\left(\cdot \right)$ is the Gamma function. On the other hand, Cole claimed that alternating current resistance and capacity measurements over a wide frequency range show that biological materials may be considered electrically equivalent to a circuit including a polarization element considered as a resistance and a capacity in series [20]. When a constant current $i$ is started through this element at time $t=0$, the potential difference across the element may be found to be:

$$e\left(t\right)=\frac{z_1{\mathit{it}}^{\alpha }}{\mathrm{\Gamma }\left(1+\alpha \right)},$$E5

where $z\left(\omega \right)=z_1{\left(\mathit{j\omega }\right)}^{-\alpha }$ is the impedance of the element, with $j^2=-1$. The phase of this element: $\phi =\mathit{\alpha \pi }/2$ allowed to provide the value of $\alpha$ for diverse biological materials. For example, it was found that $\alpha$ is ranged from 0.62 to 0.71 for frog sciatic nerve. This justifies why the power-law dynamics can occur in the membrane electrical activity. Note that Eqs. (3) and (5) are equivalent. Fractional derivatives with power-law kernel have also been used to account for the multiple timescale dynamics in neuroscience i.e. for processes which cannot be dynamically characterized with a single time constant. For example, Lundstrom et al. [21] showed that single rat neocortical pyramidal neurons adapt with a time scale that depends on the time scale of changes in stimulus statistics and that this multiple time scale adaptation is consistent with fractional-order differentiation, such that the neuron firing rate is a fractional derivative of slowly varying stimulus parameters. Subsequently, Teka et al. [22] developed a fractional-order leaky integrate-and-fire model which can reproduce upward and downward spike adaptations found experimentally by Lundstrom et al. [21]. In Ref. [22], the fractional derivation operator was applied on the membrane potential. So, the spike timing adaptation would be a dynamical manifestation of dielectric losses in the membrane. So, in order to account for the dielectric losses in the membrane, the total membrane current should be written as follows:

$$I=C_m^{\alpha }\frac{d^{\alpha }v}{{\mathit{dt}}^{\alpha }}+I_{\mathrm{Na}}+I_{\mathrm{K}}+I_l,$$E6

where $C_m^{\alpha }$ is a constant related to the capacitance of the membrane. Taking into account the ionic conductances and the gating variables, Eq. (6) is rewritten as follows:

$$C_m^{\alpha }\frac{d^{\alpha }v}{{\mathit{dt}}^{\alpha }}=I-{\overline{g}}_{\mathrm{K}}{n}^4\left(v-v_{\mathrm{K}}\right)-{\overline{g}}_{\mathrm{Na}}{m}^{3}h\left(v-v_{\mathrm{Na}}\right)-{\overline{g}}_l\left(v-v_l\right).$$E7

On the other hand, it has been shown that the dynamics of the gating variables would better be described using fractional derivatives instead of first-order derivatives [12]. Indeed, the gating dynamics is more complicated than what is traditionally assumed. Let us recall that ion channels are complex membrane proteins which provide ion-conducting, nanoscale pores in the biological membranes [23]. The gating dynamics is the spontaneous conformational dynamics of these proteins resulting in stochastic transitions between the conducting and non-conducting states also known as open and closed states of the pore, respectively. It was considered that the transitions between the open state and the closed state is a Markovian stochastic process that can be characterized by exponential residence time distributions of open and closed time intervals. However, experimental investigations had revealed that the distributions of the residence time intervals of closed states are typically not exponential. For several ion channels, these residence time distributions can be fitted by a stretched exponential function [24], or by a power-law function [25]. In Ref. [23], the time derivative of the Mittag-Leffler function which interpolates between the stretched exponential (at small time) and the asymptotic long time power-law function, was considered to show that the closed state of the ion channel can be modeled as an anomalous diffusion process over a large number of traps, described by a time-fractional diffusion equation. Accordingly, the dynamics of gating variables in the Hodgking-Huxley equation should be described by fractional-order differential equations as shown by Teka et al. [12], i.e.:

$$\frac{d^{\beta }\mathbf{z}}{{\mathit{dt}}^{\beta }}=a_{\mathbf{z}}\left(1-\mathbf{z}\right)-b_{\mathbf{z}}\mathbf{z},$$E8

where $\mathbf{z}=\left[n,m,h\right]$, and $\beta$ is a parameter related to the power-law index of anomalous diffusion. As shown in Ref. [23], the diffusion process over the substates of the closed state is a subdiffusive phenomenon, which corresponds to $0<\beta <1$.

Combining Eqs. (7) and (8), one obtains the incommensurate fractional-order Hodgking-Huxley model. Now, the variables $x$ and $y$ of the FHN model have been considered to describe the membrane voltage and the coarse-grained action of the gating variables, respectively. So, the variables $x$ and $y$ correspond to $v$ in Eq. (7) and $\mathbf{z}$ in Eq. (8), respectively. Thus, the incommensurate fractional-order Hodgking-Huxley model given by Eqs. (7) and (8) is reduced to the following incommensurate fractional-order FHN model:

$$\begin{array}{c}\epsilon D_t^{\alpha }x=y+x-\frac{1}{3}{x}^3+I,\\ D_t^{\beta }y=-x-\mathit{\delta y}+\gamma ,\end{array}$$E9

where $D_t^{\alpha }x=\frac{d^{\alpha }x}{{\mathit{dt}}^{\alpha }}$ is the fractional-order ($\alpha$-order) time derivative of $x\left(t\right)$ defined in the sense of Caputo [26] as follows:

$$D_t^{\alpha }x\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}{\int }_0^t\frac{x^{\left(1\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha }}\mathit{d\tau },$$E10

where $\alpha \in \mathrm{ℝ}$ and $0<\alpha <1$, and $\mathrm{\Gamma }\left(\cdot \right)$ is the Gamma function. The Caputo’s definition of fractional derivative is suitable for initial value problems. Note that $\underset{\alpha \to 1^{-}}{lim}{D}_t^{\alpha }x\left(t\right)=x^{\left(1\right)}\left(t\right)=\frac{\mathit{dx}\left(t\right)}{\mathit{dt}}$ [27] and then for $\alpha \to 1$ and $\beta \to 1$, the set of eqs. (9) is reduced to the classical FHN model given by Eq. (1). For $I=\delta =0,$ Eq. (9) amounts to the incommensurate fractional-order van der Pol oscillator with constant forcing studied by us in Ref. [10]. If we apply the transformation $y\to -y$ and set $\beta =\alpha$ in Eq. (9), we uncover another version of the fractional-order FHN model considered by Abdelouahab et al. [11].

3. Equivalent electrical circuit

An example of electrical circuit capable to reproduce the dynamical behavior of the fractional-order FHN model is apparently the same as the one considered by Nagmo et al. [17] as shown in Figure 1, where L is a fractional-order inductor, C is a fractional-order capacitor, and the other electronic elements are classic. Since, there is not yet a conventional manner to represent fractional-order electronic elements, we prefer to represent them as classical elements. TD is a tunnel diode whose voltage-current characteristic is given by $f\left(e\right)=i_0-\left(\left(e-e_0\right)-{\left(e-e_0\right)}^3/3K^2\right)/\rho ,$ where $i_0=f\left(e_0\right)$. The fractional capacitor and the fractional inductor are characterized by $i_C=C^{\alpha }\frac{d^{\alpha }v}{d{\tau }^{\alpha }},$ and $v_L=L^{\beta }\frac{d^{\beta }i}{d{\tau }^{\beta }},$ respectively, where $C^{\alpha }$ and $L^{\beta }$ are parameters related to their capacitance and inductance, with $0<\alpha <1$, $0<\beta <1$. Some of these coefficients can be found in Refs. [19, 28] for real capacitors and in Ref. [29] for real inductors. Applying the Kirchoff’s law, it comes out that the circuit of Figure 1 is described by the following set of fractional-order differential equations:

Figure 1.

An electrical circuit simulating the fractional-order FHN model.

$$\begin{array}{l}{C}^{\alpha }\frac{d^{\alpha }v}{d{\tau }^{\alpha }}=i+i_0-\frac{1}{\rho }\left(\left(E_0-v-e_0\right)-\frac{1}{3K^2}{\left(E_0-v-e_0\right)}^3\right)+j,\\ L^{\beta }\frac{d^{\beta }i}{d{\tau }^{\beta }}=-\mathit{Ri}-v.\end{array}$$E11

Let ${\tau }_L=L/\rho$ and ${\tau }_C=\mathit{\rho C},$ the time constants related to the dynamics of the inductor and capacitor, respectively. Let us introduce the following dimensionless variables: $t=\tau /{\tau }_L$, $x=\left(v+e_0-E_0\right)/K$, and $y=\rho \left(i+i_0\right)/K$, and use the fractional differential operator $D_t^{\alpha }$. Then, Eq. (11) can be rewritten as follows:

$$\begin{array}{c}{\mathit{\epsilon \tau }}_L^{1-\alpha }{D}_t^{\alpha }x=y+x-\frac{1}{3}{x}^3+I,\\ {\tau }_L^{1-\beta }{D}_t^{\beta }y=-x-\mathit{\delta y}+\gamma ,\end{array}$$E12

where $\epsilon ={\tau }_C/{\tau }_L={\rho }^{2}C/L$, $I=\mathit{\rho j}/K$, $\delta =R/\rho$, $\gamma =\left({\mathit{Ri}}_0+e_0-E_0\right)/K$. The difference of scales between ${\tau }_L$ and ${\tau }_C$ is at the basis of the slow-fast dynamics that results in relaxation oscillations in the FHN system behavior. Indeed, ${\tau }_C\ll {\tau }_L$, then $\epsilon ={\tau }_C/{\tau }_L\ll 1$ acts as a time scales ratio between $x$ and $y$. Without harming any generality, we will consider ${\tau }_L=1$, since the only effect of this parameter is to reinforce the relaxation that is expressed yet in $\epsilon$. Then, Eq. (12) can be rewritten without ${\tau }_L$, which amounts to the fractional-order FHN model given by Eq. (9). So, the fractional-order FHN model can be simulated with the electrical circuit in Figure 1 where the capacitor and the inductor are fractional-order electrical elements known as fractances.

4. Local stability analysis and Hopf-like bifurcation

The equilibrium points $\mathrm{E}\left(x_{\ast },y_{\ast }\right)$ of Eq. (9) are solutions of the following set of algebraic equations:

$$\begin{array}{l}{x}_{\ast }^3-3\left(1-\frac{1}{\delta }\right)x_{\ast }-3\left(I+\frac{\gamma }{\delta }\right)=0,\\ y_{\ast }=\frac{-x_{\ast }+\gamma }{\delta }.\end{array}$$E13

We will consider the case where this equation admits only one solution, i.e. for $-4{\left(1-1/\delta \right)}^3+9{\left(I+\gamma /\delta \right)}^2>0,$ according to the Cardan’s method. The fractional dynamics does not affect neither the number of equilibrium points nor their positions, but it may change their stability [26]. So, it is proper to study the stability of E in this particular context and conclude about the effect of the fractional derivatives. To do so, we will first consider the local stability of the classical integer-order FHN model. Let $\left\{\lambda \right\}$ be the eigenvalues spectrum of the Jacobian matrix ${\mathbf{J}}_{\mathrm{E}}$ of Eq. (9) evaluated at equilibrium point $\mathrm{E}\left(x_{\ast },y_{\ast }\right)$. The corresponding eigenvalues are conjugate complex numbers given by:

$${\lambda }_{1,2}=\frac{-\left(\mathit{\epsilon \delta }-1+x_{\ast }^2\right)\pm j\sqrt{\mathrm{\Delta }}}{2\epsilon },$$E14

where $j^2=-1$, and $\mathrm{\Delta }=4\epsilon \left[1-\delta \left(1-x_{\ast }^2\right)\right]-{\left(\mathit{\epsilon \delta }-1+x_{\ast }^2\right)}^2>0.$ Then, $\mathrm{E}$ is a focus, a key ingredient for the occurrence of Hopf bifurcation. According to Eq. (14), $\mathrm{E}$ is stable (for the integer-order system) if:

$$\mathit{\epsilon \delta }-1+x_{\ast }^2>0.$$E15

Thus, Hopf bifurcations occur at $x_{\ast }=x_{\mathrm{H}}^{+}=\sqrt{\mathit{\epsilon \delta }-1}$ and $x_{\ast }=x_{\mathrm{H}}^{-}=-\sqrt{\mathit{\epsilon \delta }-1}$ which correspond via Eq. (13) to two values of the stimulus current, namely $I_{\mathrm{H}}^{+}$ and $I_{\mathrm{H}}^{-}$, respectively. Figure 2 shows the bifurcation diagram computed using Matcont Matlab toolbox, around $I_{\mathrm{H}}^{+}$ for a set of parameters chosen very close to the one used by FitzHugh in its pioneering work [30], namely $\epsilon =0.1$, $\delta =0.8$, $\gamma =0.7$, that will be used all through the paper. Figure 2 shows that the Hopf bifurcation obtained for this set of parameters is subcritical.

Figure 2.

One parameter bifurcation diagram of the FHN system with respect to the stimulus current $I$, for $\epsilon =0.1$, $\delta =0.8$, $\gamma =0.7$. The solid (resp. dashed) gray line depicts stable (resp. unstable) equilibrium point; the black solid (resp. red hollow) circle markers depict the extremums of stable (resp. unstable) periodic orbits. H and LPC label Hopf bifurcation and fold bifurcation of periodic orbits also known as saddle-node bifurcation of limit cycles and limit point bifurcation for cycles, respectively.

Secondly, we consider the fractional-order FHN model with different rational orders $\alpha =m/m^{\prime }$ and $\beta =n/n^{\prime }$, with $m,n,m^{\prime },n^{\prime }\in \mathrm{\mathbb{N}}$. Let $M$ be the less common multiple of $m^{\prime }$ and $n^{\prime }$ . According to the Theorem on the stability of incommensurate fractional-order systems [26], the equilibrium point E is asymptotically stable if all the roots ${\lambda }_k$ of the following equation:

$$det\left(\mathbf{diag}\left[{\lambda }^{\mathit{M\alpha }}{\lambda }^{\mathit{M\beta }}\right]-{\mathbf{J}}_{\mathrm{E}}\right)=0,$$E16

satisfy the following condition:

$$\left|arg\left({\lambda }_k\right)\right|>\frac{\pi }{2M},\forall k,$$E17

where ${\mathbf{J}}_{\mathrm{E}}$ is the Jacobian matrix of the system evaluated at $\mathrm{E}$. The stability condition given by Eq. (17) is equivalent to the following:

$${\lambda }_k\notin \mathrm{D}=\left\{z\in \mathrm{ℂ}/z=r{\mathrm{e}}^{\mathit{j\theta }},0\le r<R\to \infty ,\left|\theta \right|\le \frac{\pi }{2M}\right\},\forall k.$$E18

And Eq. (17) is equivalent to the following characteristic equation:

$$P\left(\lambda \right)={\lambda }^{M\left(\alpha +\beta \right)}+{\mathit{\delta \lambda }}^{\mathit{M\alpha }}-\frac{1-x_{\ast }^2}{\epsilon }{\lambda }^{\mathit{M\beta }}+\frac{1-\delta \left(1-x_{\ast }^2\right)}{\epsilon }=0.$$E19

Let ${\mathrm{B}}_{\mathrm{D}}$ be the boundary of $\mathrm{D}$. According to the Cauchy’s argument Principle [31, 32], if Eq. (19) has no root on ${\mathrm{B}}_{\mathrm{D}}$ the closed curve $P\left({\mathrm{B}}_{\mathrm{D}}\right)$ encircles the origin $N$ times, where $N$ is the number of roots of Eq. (19) inside the domain $\mathrm{D}$. Accordingly, the stability condition given by Eq. (18) requires that $N=0$. Therefore, the stability condition can be resumed in the following theorem [32]:

Theorem 1: The equilibrium point E of the fractional-order FHN model is stable if $P\left({\mathrm{B}}_{\mathrm{D}}\right)$ neither encircles nor gets through the origin $\mathrm{O}$ in the complex plan.

Drawing one’s inspiration from the method of exploitation of Theorem 1 proposed in Ref. [32] and improved in Ref. [10], one can derive the following stability condition for the equilibrium point $\mathrm{E}$ of the incommensurate fractional-order FHN model:

$${\zeta }^{\alpha +\beta }cos\frac{\left(\alpha +\beta \right)\pi }{2}+{\mathit{\delta \zeta }}^{\alpha }cos\frac{\mathit{\alpha \pi }}{2}-\frac{1-x_{\ast }^2}{\epsilon }{\zeta }^{\beta }cos\frac{\mathit{\beta \pi }}{2}+\frac{1-\delta \left(1-x_{\ast }^2\right)}{\epsilon }>0,$$E20

where $\zeta$ is solution of the following equation:

$${\zeta }^{\alpha +\beta }sin\frac{\left(\alpha +\beta \right)\pi }{2}+{\mathit{\delta \zeta }}^{\alpha }sin\frac{\mathit{\alpha \pi }}{2}-\frac{1-x_{\ast }^2}{\epsilon }{\zeta }^{\beta }sin\frac{\mathit{\beta \pi }}{2}=0.$$E21

where $\zeta =r^M$, with $r\in {\mathbf{ℝ}}_{+}$. For a given value of $I$, Eq. (13) is solved and the solution $x_{\ast }$ is introduced into Eq. (21) which is solved its turn for a given set $\left(\alpha ,\beta \right)$, and the solution $\zeta$ is introduced into Eq. (20) to check the stability condition. Let us recall that in the case of the integer-order system, the stability changes via Hopf bifurcations. Now, considering its definition, a Hopf bifurcation cannot occur in a fractional-order system which cannot have exact periodic solutions on a finite time interval [33]. However, S-asymptotically T-periodic functions, can occur as solutions of a fractional-order autonomous system with fixed bounded lower terminal, instead of normal T-periodic solutions [34]. Then, the concept of Hopf-like bifurcation has been introduced to characterize the stability change of an equilibrium point giving rise to S-asymptotically T-periodic solutions [11]. As the stimulus current $I$ is the bifurcation parameter, the numerical simulation of the set of Eqs. (13), (20) and (21) shows that the stability of the equilibrium point $\mathrm{E}$ of the fractional-order FHN model switches via Hopf-like bifurcations at two critical points that we will refer to as $I_{\mathrm{H}}^{+}$ and $I_{\mathrm{H}}^{-}$, whose values depend on $\alpha$ and $\beta$ . It is worthwhile pointing out that for $\beta =\alpha \to 1$ these bifurcation points merge with those obtained above in the case of integer-order FHN model.

In the case, where the fractional dynamics appears only in the membrane potential (i.e., for $0<\alpha <1$ and $\beta \to 1$), and in the case where the fractional dynamics appears only in the gating variables (i.e., for $\alpha \to 1$ and $0<\beta <1$), the stability conditions can be derived easily from Eqs. (20) and (21). These two limiting cases are depicted in Figures 3 and 4, where we can see how the positions of the bifurcations points $I_{\mathrm{H}}^{+}$ and $I_{\mathrm{H}}^{-}$ vary with respect to the fractional derivatives orders $\alpha$ and $\beta$. Figure 5 shows the effects of the fractional-order derivatives on the stability boundaries of Figure 4(b). Figure 5 corresponds to the general case where the fractional dynamics appear in both the membrane potential and gating variables, that is, for $0<\alpha <1$ and $0<\beta <1$ . Overall, as shown in Figures 3–5, as the derivatives orders $\alpha$ and $\beta$ decrease, the region corresponding to unstable equilibrium shrinks, involving the expansion of stability regions. Thus, as expected, the fractional-order derivation enhances stability in the dynamics of the FHN system.

Figure 3.

Stability chart of the equilibrium point $\mathrm{E}$ in $\left(I,\alpha \right)$ space for $\beta \to 1$: (b) is a zoom of (a). The gray colored area is the oscillatory region, which corresponds to unstable $\mathrm{E}$. The black solid lines depict the Hopf-like bifurcations or stability boundaries $I_{\mathrm{H}}^{+}\left(\alpha \right)$ and $I_{\mathrm{H}}^{-}\left(\alpha \right)$. UE (resp. SE): Unstable (resp. stable) equilibrium point $\mathrm{E}$.

Figure 4.

Stability chart of the equilibrium point $\mathrm{E}$ in $\left(I,\beta \right)$ space for $\alpha \to 1$: (b) is a zoom of (a). The gray colored area is the oscillatory region, which corresponds to unstable $\mathrm{E}$. The black solid lines depict the Hopf-like bifurcations or stability boundaries $I_{\mathrm{H}}^{+}\left(\beta \right)$ and $I_{\mathrm{H}}^{-}\left(\beta \right)$. UE (resp. SE): Unstable (resp. stable) equilibrium point $\mathrm{E}$.

Figure 5.

Effect of $\alpha$ and $\beta$ on the stability boundary (Hopf-like bifurcation) $I_{\mathrm{H}}^{+}$.

In what comes next, keeping in mind the above local stability analysis, we will examine the oscillatory behaviors of the fractional-order FHN system.

5. Dynamical behavior and neurocomputational implications

The incommensurate fractional-order FHN model given by Eq. (9) is solved numerically thanks to the Adams-Bashforth-Moulton predictor-corrector scheme [35], with the set of parameters used for Figure 2, that is, $\epsilon =0.1$, $\delta =0.8$, $\gamma =0.7$. Unless otherwise indicated, the initial conditions are set as $x\left(0\right)=0$ and $y\left(0\right)=0$ . The fractional derivatives orders $\alpha$ and $\beta$, and the input stimulus current $I$ are assumed to be control parameters.

Let us recall that, for the chosen set of parameters, the dynamics of the classical integer-order FHN model (Eq. (1) or Eq. (9) for $\alpha \to 1$ and $\beta \to 1$) converges either to the equilibrium point (resting state) or to a limit cycle with relaxation oscillations (spiking state), depending on the strength of the stimulus $I$. Since, the transition between these two regimes occurs via a subcritical Hopf bifurcation, there is an interval of $I$ where the two attractors coexist.

In the case, where the fractional dynamics appears only in the membrane potential (i.e. for $0<\alpha <1$ and $\beta \to 1$), a narrow region of a new regime, namely phasic spiking, appears between the regions of existence of resting state and spiking state. The phasic spiking pattern is made up of a spiking phase transient to resting. The size of its region of existence in the parameter space increases with decreasing value of $\alpha$. Figure 6 shows an illustration of these three dynamical regimes. Note that the region of existence of spiking state extends beyond the stability threshold, which means that the subcritical Hopf-like bifurcation persists. Keeping the value of $\alpha \in \left]0,1\right[$, the value of $\beta$ is reduced a bit, say from 1 to 0.9. The bifurcation scenario changes a lot. When the system loses stability, MMOs take place (see Figure 7) and when the value of $I$ decreases further, there is a transition to spiking state. For high values of $\alpha$, at the transition between resting state and MMOs, a narrow window of small amplitude oscillations (see Figure 8) appears. A very narrow region of another type of phasic spiking appears, for relatively small values of $\alpha$, at the transition from resting state to MMOs. This new type of phasic spiking is made up of transient MMOs to resting (see Figure 8). The lifetime and the number of spikes of this transient MMOs regime increases as the strength of the stimulus current goes away from the bifurcation point.

Figure 6.

Time series $-x\left(t\right)$ showing the dynamical regimes exhibited by the FHN model with fractional-order dynamics only in the membrane potential (i.e. $\beta \to 1$) for $\alpha =0.8$, and different values of the stimulus $I$: (a) $I=-0.400$: Spiking; (b) $I=-0.345$: Phasic spiking; (c) $I=-0.330$: Resting.

Figure 7.

Time series $-x\left(t\right)$ showing some MMOs patterns exhibited by the fractional-order FHN model for $\alpha =0.8$, $\beta =0.9$, and different values of $I$: (a) $1^{15}$ pattern for $I=-0.4146$; (b) $1^3{1}^2$ pattern for $I=-0.4290$; (c) $1^2$ pattern for $I=-0.4320$; (d) $1^1{2}^1$ pattern for $I=-0.4472$.

Figure 8.

Time series $-x\left(t\right)$ showing some dynamical regimes exhibited by the fractional-order FHN model: (a) small amplitude oscillations for $\alpha =0.985$, $\beta =0.9$, and $I=-0.3824$; (b) phasic spiking for $\alpha =0.6$, $\beta =0.9$, and $I=-0.427$.

Regular MMOs have been often referred to as $L^s$ patterns, where $L$ and $s$ are the number of large amplitudes and the number small amplitudes in one pattern, respectively. For example, Figure 7 shows MMOs with $1^{15}$ pattern in (a), and $1^2$ pattern in (c). When the value of the control parameter $I$ varies, MMOs with complex patterns develop at the transition between two MMOs states with regular patterns; for example, Figure 7 shows MMOs with $1^3{1}^2$ pattern in (b); and $1^1{2}^1$ pattern in (d). The $1^3{1}^2$ MMOs develops between $1^3$ MMOs and $1^2$ MMOs, and $1^1{2}^1$ MMOs develops between $1^1$ MMOs and $2^1$ MMOs. Figure 7 also shows that when the value of $I$ is close to the bifurcation point $I_{\mathrm{H}}^{+}\approx -0.41185$, we have $1^s$ MMOs – the value of $s$ decreases for decreasing value of $I$. As the value of $I$ decreases further and the value of $s$ reaches 1, one obtains complex MMOs patterns with many large amplitude oscillations, which finally leads to $L^1$ MMOs (where $L$ increases with decreasing value of $I$). The MMOs patterns shown in Figure 7 are identical to classical folded nodes-induced MMOs observed in integer-order systems, also known as canard generated MMOs [16].

In order to illustrate the bifurcation scenarios described above, we map all the dynamical regimes in the parameter space as shown in Figure 9. As mentioned above, the Hopf-like bifurcation obtained there is subcritical as shown in Figure 9(a). Indeed, there is a parameter region where oscillatory states coexist with resting state. However, Figure 9(b) show that when the value of $\beta$ decreases, the limit between oscillatory states and resting state is marked clearly by the stability threshold obtained by local stability analysis.

Figure 9.

Dynamical regimes maps in $\left(I,\alpha \right)$ space superimposed on the stability boundaries $I_{\mathrm{H}}^{+}\left(\alpha \right)$ obtained analytically: (a) $\beta \to 1$; (b) $\beta =0.9$. Gray colored region: Spiking; cyan colored region: MMOs; magenta colored region: Small amplitude oscillations (SAOs); red colored region: Phasic spiking (PS), white region: Resting. The square markers in (a), circle markers and diamond marker in (b) show the parameters used for Figures 6–8, respectively.

We now examine the case where the fractional dynamics appears only in the gating variables, that is, for $\alpha \to 1$ and $0<\beta <1$ . The resting state exists for high values of $I$. As the value $I$ decreases the equilibrium state loses stability and MMOs take place. At the transition between resting state and MMOs, a very narrow window of small amplitude oscillations appears. As the value of $I$ decreases further, there is a transition from MMOs to spiking state. Keeping $0<\beta <1$, the value of $\alpha$ is reduced a bit, say from 1 to 0.9. The bifurcation scenario does not change significantly. Finally, we map all the dynamical behaviors in the parameter space as shown in Figure 10. One can notice that the result of the global dynamics analysis obtained numerically agrees very well with the local stability analysis result. This figure shows that the domain of existence of MMOs widens as the value of $\beta$ decreases. The major difference between the cases depicted in the subsets (a) and (b) of Figure 10 is that the region of existence of small amplitude oscillations shrinks with decreasing value of $\alpha$.

Figure 10.

Dynamical regimes maps in $\left(I,\beta \right)$ space superimposed on the stability boundaries $I_{\mathrm{H}}^{+}\left(\beta \right)$ obtained analytically: (a) $\alpha \to 1$; (b) $\alpha =0.9$. Gray colored region: Spiking; cyan colored region: MMOs; magenta colored region: Small amplitude oscillations (SAOs); white region: Resting. The circle marker in (b) show the parameters used for Figure 11.

These phase diagrams in Figures 9 and 10 show the rich variety of dynamical behaviors of the fractional-order FHN neuron model.

For low values of $\beta$ (say, lesser than about 0.7) one observes another type of MMOs for which small oscillations start with very low amplitude which then grows slowly before the oscillations enter the spiking phase (see Figure 11). This last type of MMOs is identical to singular Hopf bifurcation-induced MMOs observed in integer-order systems [16]. Note that unlike the integer-order FHN model for which in the unstable equilibrium region, the dynamics of the system converges to a unique limit cycle no matter the initial conditions, the dynamics of its fractional-order counterparts is sensitive to initial conditions [11]. The same phenomenon was observed in the fractional-order van der Pol oscillator with constant forcing [10] which is closely related to the FHN model. So, changing the initial conditions, it is possible to uncover new dynamical states. For example, we have change the set of initial conditions $\left(x\left(0\right),y\left(0\right)\right)$ in Figure 11 from $\left(0,0\right)$ to $\left(0,1\right)$, and we have found that the active state is preceded by a static-like regime. This static-like transient state was found for the first time in the context of dynamical systems by us, while studying a fractional-order van der Pol oscillator with constant forcing [10]. In the context of neuroscience, this quiescent transient regime is known as first spike latency. Teka et al. [22] found that the fractional derivation also induces the occurrence of first spike latency in a leaky integrate-and-fire model, and that this first spike latency is reinforced by decreasing value of the fractional derivative. On the other hand, still regarding the type MMOs shown in Figure 11, the instantaneous spike frequency (the inverse of the duration of the corresponding interspike interval), also known as the firing rate, increases as a function of the interspike interval number and approaches an asymptotic value after a certain number of intervals. Also, the increasing rate of the firing rate decreases with decreasing value of the fractional derivatives orders. This property of neurons known as upward spike timing adaptation, observed experimentally [21], has proved to appear only in fractional-order models [21, 22].

Figure 11.

Time series $-x\left(t\right)$ showing two different behaviors of the fractional-order FHN model for $\alpha =0.9$, $\beta =0.6$, and $I=-0.540$ and different initial conditions: (a) $\left(x\left(0\right),y\left(0\right)\right)=\left(0,0\right)$; (b) $\left(x\left(0\right),y\left(0\right)\right)=\left(0,1\right)$.

It was found in Ref. [10] that the lifetime of this static-like transient and the pseudo-period (corresponding to the interspike interval) of the asymptotic MMOs increase exponentially with the closeness to the Hopf bifurcation point. The same result is found in the present work. Figure 12 shows the behavior of the first spike latency (FSL) and of the fiftieth interspike interval (ISI) – the period of the asymptotic MMOs after spike timing adaptation - when the value of the stimulus current $I$ varies, coming close to the bifurcation point $I_{\mathrm{H}}^{+}$. The results from the direct numerical simulation of Eq. (9) (depicted in Figure 12) have been fitted using the following exponential functions in EzyFit Matlab toolbox:

Figure 12.

First spike latency (FSL) and interspike interval (ISI) versus stimulus current $I$, for $\alpha =0.9$ and $\beta =0.6$ . With $I_H^{+}\approx 0.504497$.

$$\mathit{FSL}\left(I\right)=aexp\left(\frac{b}{I-I_{\mathrm{H}}^{+}}\right)+c,\mathit{ISI}\left(I\right)=a^{\prime }exp\left(\frac{b^{\prime }}{I-I_{\mathrm{H}}^{+}}\right)+c^{\prime }.$$E22

For the two cases of fitting, the value of the correlation coefficient is over 0.999, indicating good fittings. Thus, this confirms the exponential growth of the first spike latency and of the interspike interval. Furthermore, as shown in Figure 13, the features of the first spike latency and of the interspike interval are reinforced by decreasing values of the orders of fractional derivatives.

Figure 13.

First spike latency (FSL) versus displacement of the stimulus current $I$ from the Hopf-like bifurcation point $I_{\mathrm{H}}^{+}$, for $\alpha =0.9$ and different values of $\beta$. Note that $I_{\mathrm{H}}^{+}\left(\beta =0.6\right)\approx -0.504497$ and $I_{\mathrm{H}}^{+}\left(\beta =0.8\right)\approx -0.431551$.

Although the specific role of MMOs among the plethora of slow-fast dynamical behaviors occurring in neural systems has not yet been determined, it has been nevertheless suggested that the dynamical richness of a neuron expands its computational capacity by increasing its coding capacity [12, 15, 36]. On the other hand, it has been suggested that first spike latency could code for stimulus recognition in several sensory systems [37, 38, 39, 40]. Indeed, the variation of first spike latency with stimulus parameters contains considerable information about those parameters [40]. The first spike latency has also been suggested as a source of information for accurate decisions [22]. In addition, the behavior of the firing rate (or interspike interval) of neurons also provides information on the stimulus statistics [21]. The fractional-order FHN model is a mathematically simple model with complex dynamics features, thus increasing the amount of information that can describe the input. So, the fractional FHN model is more computationally efficient than its integer-order counterpart.

6. Conclusion

In this work, the dynamical behavior of an incommensurate fractional-order FitzHugh-Nagumo model of neuron has been investigated in details. First of all, the considered fractional-order model has been derived from the fractional-order Hodgkin-Huxley model obtained taking advantage of the powerfulness of fractional derivatives in modeling the dielectric losses in cell membranes, and the anomalous diffusion of particles in ion channels. Then, it has been shown that the fractional-order FitzHugh-Nagumo model can be simulated by a simple electrical circuit. Then, the local stability of the incommensurate fractional-order model has been studied with a particular attention granted to the effect of the fractional derivatives. It has been found that the fractional derivatives enhance the local stability of the model. At last, the dynamical behavior of the fractional-order models has been explored numerically, which has confirmed the results of local stability. It has been found that the fractional-order FitzHugh-Nagumo exhibits a lot of complex dynamical features that are not observed in the behavior of its integer-order counterpart, and that cannot be observed in integer-order two-variable systems. Among others, the fractional-order model exhibits mixed mode oscillations, phasic spiking, first spike latency, and spike timing adaptation. These complex features of the dynamical behavior of the fractional-order model increase the computational capacity of the FitzHugh-Nagumo model. An outlook of this work would be the study of the mechanism(s) underlying the formation of mixed mode oscillations and first spike latency as effects of fractional derivation in two-variable systems.

Conflict of interest

The author declares no conflict of interest.