A fractional-order improved FitzHugh–Nagumo neuron model

Understanding the underlying mechanics of information processing and cognitive function is made easier by knowledge of the biological dynamics of neurons. In neuron-based application engineering, a lot of electrical activity in response to an externally supplied stimulus is crucial. The Hodgkin–Huxley neuron model,^[1,2] the Morris–Lecar neuron model,^[3,4] the Chay neuron model,^[5,6] the Hindmarsh–Rose neuron model,^[7,8] and the FitzHugh–Nagumo (FHN) neuron model^[9,10] are just a few of the many neuron models that have been constructed to date to mimic various biological neuron dynamics. Generalized from the Hodgkin–Huxley model, the FHN model has good theoretical and numerical simulation efficiency for neuron electrical activity.

There are several studies that have been reported on the various types of neuron models. In Ref. [11], the authors proposed a nonlinear fitting algorithm for analogue multiplierless application of the Hindmarsh–Rose neuron model. In Ref. [12], the authors analyzed a digital multiplierless realization of a two-coupled biological Hindmarsh–Rose neuron model. In Ref. [13], Cai et al. proposed a digital multiplierless implementation for a nullcline-based piecewise linear Hindmarsh–Rose neuron model. In Ref. [14], the authors studied a nullcline-based control technique for PWL-shaped oscillators. In Ref. [15], the authors analyzed the types of bifurcations of FHN maps. In Ref. [16], a study on the dynamics of the FHN neuron model defined with Levy noise and Gaussian white noise was proposed. In Ref. [17], a study on the weak periodic signal detection by sine-Wiener-noise-induced resonance in the FHN neuron was proposed. In Ref. [18], the authors proposed a modified FHN neuron model with multiplier-free implementations. In Ref. [19], the authors studied temperature-induced logical resonance in the Hodgkin–Huxley neuron. In Ref. [20], Guo et al. studied Hopf bifurcation and phase synchronization in memristor-coupled Hindmarsh–Rose and FHN neurons with two time delays.

Fractional calculus is one of the highly used tools of applied mathematics while dealing with the modeling of real-life problems.^[21–23] Various types of fractional-order operators are available, of which Riemann–Liouville (RL) and Caputo are the most commonly used. Recently, a fractional derivative with an exponential delay-type kernel was proposed in Ref. [24]. In Ref. [25], Atangana and Baleanu proposed a derivative with a Mittag–Leffler kernel. There are several types of Caputo derivatives that have been derived in the recent studies.^[26–28] In Ref. [29], the authors proposed a generalized form of the Caputo fractional derivative. Some recent applications of the generalized Caputo fractional derivative have been seen in epidemiology,^[30,31] ecology,^[32,33] psychology,^[34] and neuron dynamics.^[35]

Solving the differential equations in terms of fractional derivatives is one of the research interests of many researchers. Several numerical methods have been recently proposed to derive the numerical solutions of different types of fractional differential equations (FDEs). In Ref. [36], the authors proposed a predictor–corrector (PC) scheme to solve Caputo FDEs. In Ref. [37], the L1-based predictor–corrector method was proposed to solve Caputo-type problems. Later, this method was used to solve delay-type problems in Ref. [38]. In Ref. [39], the authors proposed a high-order PC scheme on graded meshes to solve FDEs. In Ref. [29], the researchers proposed a PC method for generalized Caputo-type problems. A PC scheme for generalized Caputo problems with time delay was derived in Ref. [40]. The authors in Ref. [41] proposed a neural network-based scheme to solve generalized Caputo-type problems. An efficient finite-difference PC method was introduced in Ref. [42]. Recently, in Ref. [43], the authors proposed an L1-PC method for solving generalized Caputo FDEs.

In the literature, fractional derivatives have been significantly incorporated into neuron models to simulate complex information processing and capture memory effects. The fractional-order models offer additional design flexibility by providing adjustable fractional orders compared to integer-order models. Sacu in Ref. [44] proposed a fractional-order FHN neuronal model with experimental analysis. In Ref. [45], the authors proposed a fractional-order Hindmarsh–Rose neuronal model. In Ref. [46], the authors analyzed two- and three-dimensional fractional-order Hindmarsh–Rose neuron models. In Ref. [47], the authors proposed some hidden dynamics in a fractional-order memristive Hindmarsh–Rose model. Mondal et al.^[48] studied the firing patterns of a fractional-order FitzHugh–Rinzel bursting neuron model. In Ref. [49], the authors analyzed spiking and bursting patterns of the fractional-order Izhikevich model. In Ref. [50], the authors studied FPAA-based implementation of fractional-order chaotic oscillators using first-order active filter blocks. In Ref. [51], Malik and Mir proposed a discrete multiplierless implementation of the fractional order Hindmarsh–Rose model. In Ref. [52], Tolba et al. performed the analyses of synchronization and FPGA realization of the fractional-order Izhikevich neuron model. In Ref. [53], the authors studied low-voltage low-power integrable CMOS circuit application of integer- and fractional-order FHN neuron models. In Ref. [54], the authors studied the spiking and bursting of a fractional modified FHN neuron model. In Ref. [55], dynamical modeling and equilibrium stability of a fractional discrete biophysical neuron model were provided. In Ref. [35], Kumar et al. proposed a generalized Caputo-type fractional-order neuron model under the electromagnetic field.

In this paper, we propose a fractional version of an improved FHN neuron model^[18] in the sense of generalized Caputo fractional derivatives. The motivation behind proposing a fractional-order improved FHN neuron model is that such a system can provide a more nuanced description of the process with better understanding and simulation of the neuronal responses by incorporating memory effects and non-local dynamics, which are inherent to many biological systems. Also, such a system has limitless memory, which can improve the system’s controllability for a wide range of real-world problems with significant applications. The reason for preferring the generalized Caputo fractional derivative compared to other aforementioned derivatives is that this derivative has extra parameters involved with the order of the operator, which offers additional design flexibility by providing two adjustable fractional orders compared to integer-order and other fractional-order operators. Later, the proposed model is solved by utilizing a recently proposed L1 predictor–corrector scheme given in Ref. [43]. The proposed fractional generalization of the integer-order FHN neuron model^[18] along with a recently proposed numerical methodology is the main contribution of this work.

This paper is designed as follows. In Section 2, we provide some definitions and the DGJ algorithm. Section 3 provides the novel fractional-order FHN neuron model. In Section 4, the necessary qualitative and quantitative analyses of the proposed model are given, where we prove the existence of a unique solution and then drive the numerical solution of the model. The graphical simulations are given in Section 5. We conclude our results in Section 6.

Definition 1^[56] The generalized Caputo fractional integral of an absolutely continuous function f:[t₀,T]→ℝ of order η is given by

Definition 2^[29,56] The generalized Caputo fractional derivative of an absolutely continuous function f:[t₀,T]→ℝ of order η is given by

For 0 < η≤1,

Consider the equation

containing a known function Λ and a nonlinear operator N(f).

Let the solution of Eq. (3) be written by

Then, N(f) is decomposed as

So, we get

Comparing both sides, we obtain

Therefore,

The truncation of the equation to m finite number of terms is taken by

In Ref. [18], the authors proposed an improved FHN neuron model in a dimensionless mathematical form by incorporating a nonlinear function with an N-shaped curve and multiplier-free incorporation. The model is defined as follows:

where x and y are the membrane potential and sodium gating variable, respectively; parameters a and b are constants with respect to the equilibrium states of potassium potential and sodium potential, respectively; and c is a time constant for controlling the sodium gating rate. I ≈ I_m sin(2πft) is the externally applied stimulus to the neuron, where f is the frequency, I_m is the amplitude, and t is system time.

As we know, classical systems do not contain memory effects because of their local nature. Fractional-order operators can be incorporated into a neuron model to capture the intricate memory and information processing functions of the human brain. Therefore, the fractional-order version of the above-given FHN neuron model (4) is given as follows:

where t0CDtη,ρ is the generalized Caputo fractional derivative of order 0 < η ≤ 1 with an extra variable 0 < ρ ≤ 1 with respect to time t. The significance of other variables and parameters is the same as in the integer-order model (4).

We first prove the existence of a unique solution with the help of some well-known results. Then, we will derive the numerical solution of the given FHN model using a recently proposed L1-PC method.

Here we provide the existence of a unique solution for the following initial value problem (IVP) which represents the above given model (5):

where 0 < η≤ 1, 0 < ρ ≤ 1, f = (x, y), f₀ = (x₀, y₀), and F = (2tanh(x)–x–y + I_msin(2πft), c(x+a–by)).

The Volterra integral equation for the IVP (6) is given by^[13]

Theorem 1 (existence)^[30] For 0 < η ≤ 1, f₀ ∈ ℝ, T^* > 0, J > 0, let us consider the set S:={(t,f):t∈[0,T*],|f−f0|≤J} and a continuous function F:S→ℝ. Let M:=sup(t,f)∈S|F(t,f)| and

Then, there exists a function f∈C[0,T] which satisfies the IVP (6).

Theorem 2 (uniqueness)^[30] If for the given parameters and the set S in Theorem 1, the continuous function F:S→ℝ satisfies the Lipschitz constraint for variable f, i.e.,

where V > 0 is a constant independent to t, f₁, and f₂. Then, there exists a unique solution f∈C[0,T] for the IVP (6).

Here we propose the L1-PC method^[43] to derive the numerical solution of the proposed FHN neuron model (5). The method is based on the L1-based discretization algorithm and the spline interpolation scheme. Let us consider the initial value problem

where 0 < η < 1, 0 < ρ ≤ 1, f = (x, y), f₀ = (x₀, y₀), and F = (2tanh(x)–x–y+I_msin(2πft), c(x+a–by)).

The interval [t₀,T] is discretized into ϱ parts with mesh points

where h=Tρ−t0ρϱ+1. From the substitution z = s^ρ, discretization of the interval, and applying the forward discretization to the derivative, we get

Then, the function F1(z1ρ,f(z1ρ))=(z1ρ)1−ρ is interpolated using linear splines with the nodes pointed at tνρ. Utilizing the trapezoidal quadrature rule with the weight function (t_ϱ + 1 – ⋅)^−η, we get the expression

Using the integral (10) in Eq. (11), we get

The coefficients a_ϱ,ν and b_ϱ,ν are defined by

Using a_ϱ,ν and b_ϱ,ν in Eq. (12) and rearranging the terms, we get the algebraic expression

Now according to the DGJ algorithm, we have

By applying the DGJ algorithm, the approximations for f(t_ϱ+1) = f_ϱ+1 are derived by

Therefore, f(t) at t = t_ϱ+1 is defined by

Finally, the expression (14) gives us the following L1-PC algorithm:

where fϱ+1P and zϱ+1P are the predictors and fϱ+1C is the corrector term.

The algorithm (15) can be written in a compact form as follows:

where C_ϱ,ν is given by

From the Eq. (13), it is clear that a_ϱ,ν and b_ϱ,ν are positive. To verify the positivity of C_ϱ,ν, we have

Also, if tϱ+11−ρ>tϱ1−ρ>tϱ−11−ρ. Thus, for ν = 1,…,ϱ, we get

So, C_ϱ,ν is positive. Also, C_ϱ,0 and μ are positive because a_ϱ,ν and b_ϱ,ν are positive. Thus, C_ϱ,ν/μ > 0.

Therefore, the numerical solution of the FHN neuron model (5) is derived as follows:

Here we provide the error estimation for the given scheme.

Theorem 3^[43] Let the generalized Caputo derivative of f(t) at t = t_ϱ+1 be defined by t0CDtη,ρf(tϱ+1). Then,

where C is a positive constant depends only on η and ρ such that

Remark 1 For ν = 1,…,ϱ, we have

Similarly, for ν = 0, we obtain

Also,

Using the bounds for C_ϱ,ν and μ given in Eqs. (19) and (20), we get

where M = C_ϱ/η is a positive constant.

So, C_ϱ,ν/μ > 0 is bounded and ∑ν=0ϱCϱ,ν/μ=1.

Theorem 4 Let f(t_n) be the exact solution, fnP be the approximate solution derived from the given method, and l=maxn,ν{μCn,ν}. Then,

where C_p = 2Cl is a positive constant.

Proof Here, we use mathematical induction. Let enP=f(tn)−fnP. From Eqs. (16) and (18), we have

Similarly, we prove that the expression is true for n.

where C_p = 2Cl is a positive constant.

Remark 2 Using Eq. (21), it can be defined that |enP|≤2CTη−1h3−η/M, where M = η/max{η,6–4η} ≤ 1 and C is a positive constant given in Eq. (18). Also, it converges to zero when h →∞.

Theorem 5 Let f(t) be the exact solution of the given IVP (9) and F(t, f) follow the Lipschitz constraint with respect to f with a Lipschitz constant L. Then,

where fnC is the approximate solution derived from the given method at t = t_n, for 0 < η < 1, and C_c is a positive constant.

Proof Consider |enC|=|f(tn)−fnC|. We have

Then, we obtain enC as follows:

where C_c is a positive constant.

Theorem 6 Let fnC and gnC represent the numerical solutions of the IVP (9) obtained by the proposed algorithm (15) or (16). Let F(t, f) follow the Lipschitz constraint with respect to f with a Lipschitz constant L, then the given method is stable.

Proof From our method

Now suppose F(t,f) follows Lipchitz constraint with respect to f with Lipchitz constant L, then Eq. (23) is expressed by

Considering only first term of Eq. (24), we get

Using the Gronwall Inequality^[38] and the bound given in Eq. (21) in the expression (25), we get

where H is a constant. Using Eq. (26) in Eq. (24), we get

From the Lipchitz condition for F, we have

Finally, using Eq. (28) in Eq. (27), we have

where

is a constant.

Hence, the required result is obtained.

Now we explore the dynamics of the proposed FHN neuron model (4) by performing graphical simulations with the help of the numerical algorithm (17). The proposed numerical solution (17) is coded in Mathematica 13.2. The model parameters a, b, and c are fixed as a = 0.70, b = 0.80, c = 0.10. In Fig. 1, the bifurcation diagram for frequency f ∈ [0.01, 0.30] vs. membrane potential x (Fig. 1(a)) and Lyapunov exponents (Fig. 1(b)) at I_m = 1, η = 0.98, and ρ = 0.99 are plotted. We notice that periodic spiking patterns are generated in the range f ∈ (0, 0.1). The chaotic patterns exist three times in a narrow parameter range of f. The quasi-periodic patterns are triggered nearly at five points of the frequency f. In Fig. 2, the bifurcation diagram for amplitude I_m ∈ [1, 2] vs. membrane potential x (Fig. 2(a)) and Lyapunov exponents (Fig. 2(b)) at f = 0.16, η = 0.98, and ρ = 0.99 are plotted. In this case, we obtain only periodic spiking patterns. In Fig. 3, the bifurcation diagram for fractional-order η ∈ [0.95, 1] vs. membrane potential x at f = 0.16, I_m = 1, and ρ = 0.99 is plotted.

In Fig. 4, we plot the dynamics of the system for f = 0.01 and I_m = 1 at fractional order η = 0.98 and parameter ρ = 0.99. Here, we achieve periodic spiking electrical activities. The same periodic spiking electrical patterns are captured in Fig. 5 for f = 0.01 and I_m = 2. In Fig. 6, for f = 0.171 and I_m = 1, we obtain the chaotic spiking electrical patterns. In Fig. 7, for f = 0.23 and I_m = 1.6, the quasi-periodic electrical activity is obtained.

We performed further graphical simulations to check the influence of fractional order on the given system. In Fig. 8, we plotted the dynamics of the system for f = 0.01, I_m = 2, and ρ = 0.99 at two different fractional order values, η = 0.99 and η = 0.96. Here, we see that periodic spiking electrical patterns exist for both fractional order values.

Similarly, in Fig. 9, we plotted the dynamics of the system for f = 0.23, I_m = 1.6, and ρ = 0.99 at two different fractional order values, η = 0.99 and η = 0.96. Here we achieve quasi-periodic electrical patterns for both fractional order values.

From the given graphical simulations, we have seen that the proposed system shows periodic, quasi-periodic, and chaotic electrical patterns at various values of the amplitude, frequency, and fractional order. The generalized Caputo fractional derivative offered a significant advantage by providing a broad representation of the memory and dynamic behavior of the proposed FHN neuron model. In comparison, the integer-order model often fails to capture the long-range temporal dependencies and the complex, history-dependent nature of neuronal activity. Therefore, the given fractional FHN neuron model provided a more nuanced description of the process with better understanding and simulation of the neuronal responses by incorporating memory effects and non-local dynamics, which are inherent to many biological systems.

In this article, we have proposed a fractional-order improved FHN neuron model in terms of a generalized Caputo fractional derivative. We have discussed the existence and uniqueness analyses and have derived the numerical solution of the model using a recently proposed L1-predictor–corrector method. We have provided the error and stability analyses of the given method. We have performed several graphical simulations to justify that the proposed FHN neuron model generates periodic spiking patterns, chaotic patterns, and quasi-periodic patterns. The role of the fractional order has been explored. The proposed fractional FHN neuron model provided a more nuanced description of the process with better understanding and simulation of the neuronal responses by incorporating memory effects and non-local dynamics, which are inherent to many biological systems. The reason behind using the generalized Caputo derivative over the standard Caputo derivative was that it offered greater flexibility and accuracy in the model with memory effects because it contained an extra parameter ρ along with the fractional order η. While the Caputo derivative captures memory by considering a fixed fractional order, the generalized Caputo derivative extends this by allowing the order of differentiation to vary with two parameters. In the future, the given model can be used for a real-time circuit experiment.