Analytical and numerical investigation of the Hindmarsh-Rose model neuronal activity

1.   Introduction

Since mathematical models are used to recreate some observed real-world issues, mathematical biology has found utility in a variety of fields. Numerous times, the dynamical mechanisms demonstrated by the transmission of infectious diseases have been faithfully simulated by mathematical models. Numerous authors have looked into the unique dynamic biological processes that occur in neurons. There have been some proposed mathematical models ^[1,2]. These models are used to illustrate the characteristics of particular nervous system cells that can produce piercing electrical currents transverse to their cell membrane ^{[3,4,5,6,7,8,9]}. In the literature, action potentials, also known as spikes, are noted for having a duration of roughly one millisecond. According to the material that is currently accessible, spikes are transferred from the conveyance neuron to numerous other neurons via axons and synapses. Spiking neurons are important in this process since they are the units that process information for coordinated nervous system function. Different types of mathematical representations of spiking neuron processes exist. For instance, Hodgkin-Huxley models are utilized to explain how the input current and ion channel activation affect the membrane voltage ^[10,11]. Integrate-and-fire models, which explain the membrane voltage as a function of input current and forecast the spike times without detailing the biophysical mechanisms that affect an action potential's time course, are simpler mathematically. Researchers have used both deterministic and stochastic approaches in the pursuit of accurately describing these processes. In addition to these ideas, differential operators have been extensively employed to incorporate other nonlocal characteristics into mathematical models. The Hindmarsh-Rose model, for instance, was developed to study the spiking-bursting behaviors of the membrane prospective that were noticed in experiments with a single neuron. These models have been presenting some chaotic phenomena. This model describes the movement of ions across the membrane via ion channels as a set of nonlinear differential equations with dimensionless unknowns. We want to perform various analyses to look into this model more thoroughly.

1.1.   Definitions and important preliminaries

In the literature, we know that the first definition of a differential operator was suggested by Leibniz and Newton. If a function $f\in C[a, b]$, they suggested that the rate of change of the function can be obtained with

Then we can obtain the derivative of the function of $f.$ The above formula gave birth to differential calculus and its inverse operator known as integral calculus. This operator was developed by many researchers. Now, we present the definitions for fractional differential and integral operators with local and nonlocal kernels introduced in ^[12,13,14].

First, we give the definition of the Riemann-Liouville fractional derivative of the function $f(t)\in C(0, T)$ as

with $\alpha$ belonging to $0 < \alpha\leq1.$

The definition of the Caputo fractional derivative of the function $f(t)\in H^{1}\left[0, T\right]$ is given by

with $\alpha$ belonging to $0 < \alpha\leq1.$

The Caputo-Fabrizio fractional derivative of the function $f(t)\in H^{1}\left(0, T\right)$ is given by

with $\alpha$ belonging to $0\leq\alpha\leq1.$

The Atangana-Baleanu fractional derivative in the Riemann-Liouville sense of the function $f(t)\in H^{1}\left(0, T\right)$ is given by

where $AB(0) = AB(1) = 1$ and $AB(\alpha) = 1-\alpha+\frac{\alpha}{\Gamma(\alpha)}$ with $\alpha$ belonging to $0\leq\alpha\leq1.$

Finally, the Atangana-Baleanu fractional derivative in the Caputo sense of the function $f(t)\in H^{1}\left(0, T\right)$ is given by

where $AB(0) = AB(1) = 1$ and $AB(\alpha) = 1-\alpha+\frac{\alpha}{\Gamma(\alpha)}$ with $\alpha$ belonging to $0\leq\alpha\leq1.$

The associated integral operators of the above fractional derivatives are defined by

2.   Analysis of Hindmarsh-Rose model of neuronal activity

To accommodate readers that are not acquainted with the Hindmarsh-Rose model, we will recall that the model was designed to investigate the spiking-bursting conduct of the membrane potential perceived in experimentations completed with a single neuron. To construct this model, some dimensionless units were considered $v(t)$, where this variable or function is the membrane potential. Two additional variables were considered, including $u(t)$ and $w(t)$. They take into account the transport of ions through the membrane via the ion channels. In particular, the function $u(t)$ measures the transport of sodium and potassium ions through fast ion channels. The function $w(t)$ represents the adaptation current. The mathematical model replicating such a dynamic is given below ^[1]:

with the initial conditions

The final equation, however, is crucial because it accounts for a wide range of dynamic membrane potential processes that are repeated by the variable. More crucially, it takes into consideration unpredictable patterns, which inevitably result in a chaotic attractor. This innate characteristic allows the model, despite its simplicity, to give a solid qualitative description of a number of processes that are present in a variety of real-world circumstances.

2.1.   Equilibrium points of system

In this subsection a detailed analysis of equilibrium points are presented. To obtain the equilibrium points of the system, we must solve the following system:

From the second equality we will get

Also from the third equation we will get

Then we will put two of them into the first equation with dependence on the variable $v^{\ast}$; we will get

For instance, by choosing $a = 1$, $b = 3,$ $c = 1,$ $d = 5$, $I = 2,$ $s = 4$ and $v_{rest} = -1.6$ and solving the cubic equation, we will get three equilibrium points given as

2.2.   Determination of $C_{0}$

In this section, we will discuss the chaotic number of the Hindmarsh-Rose model while considering the next-generation matrix method. Our system is given as

Let us divide the system into two parts. $f$ is associated with the nonlinear part of the system and $v$ is associated with the linear part of the system as below:

So we have the following matrices:

and

From the above matrices we will calculate $F$ and $V$ which are partial derivatives of $f$ and $v.$

The matrices $K = F.V^{-1}$ is called a next-generation matrix of the system. $C_{0}$ (chaotic number) is the spectral radius of the next-generation matrix of $K.$ We noted that $\rho(K)$ is the eigenvalue with maximum modules in the spectrum of $K.$

then we have

Let us calculate the eigenvalues from the equality $det(F.V^{-1}-\lambda I) = 0$ then, we have

we have three eigenvalues:

for the positive eigenvalue of $\lambda_{3}$ we must consider the following condition,

So the chaotic number for the system is given as

So ıf we consider $b > d+\frac{3av}{2}$ then we get chaos.

2.3.   Global stability of the equilibrium point via Lyapunov function

Theorem 1: If (chaotic number) $C_{0}\geq1$, then the equilibrium point $E^{\ast}(v^{\ast}, u^{\ast}, w^{\ast})$ is globally asymptotically stable.

Proof: We prove this using the idea of the Lyapunov function. We start by defining the Lyapunov function associated with the system as below:

By taking the derivative of the Lyapunov function with respect to $t$, we get

Replacing each rate of change with respect to time for each class by it expression, we obtain the following equation

After multiplying all items with each other and dividing the last equality's negative and positive parts then we write

here

and

Therefore if

3.   Existence and uniqueness analysis of system

In this section, we present some results about the existence and uniqueness of the system equations describing the neuronal activity. To show this, we define the norm

Here, we consider Banach space. Now, we give the following theorem, which is given for verifying the linear growth and Lipschitz condition properties ^[15].

Theorem 2: Assume that there exist six positive constants $l_{i}$ and $\overline{l}_{i}$ such that

ⅰ) $\forall i\in\{1, 2, 3\},$

ⅱ) $\forall i\in\{1, 2, 3\},$ $\forall(t, x)\in R^{3}\times\lbrack0, T],$

We now consider model as below;

First, we start with the function of $F_{1}(t, v(t)).$ Then, we will show that

Here, we remember the following norm:

then, we have $v, v_{1}\in R^{2}$, $t\in\lbrack0, T]$ and

After the above step, we also assume that $\forall t\in\lbrack0, T]$ and that there exist three positive constant $M_{1}, M_{2}, M_{3} < \infty$ such that $\left\Vert v\right\Vert _{\infty} < M_{1},$ $\left\Vert u\right\Vert _{\infty} < M_{2}$ and $\left\Vert w\right\Vert _{\infty} < M_{3}.$ Now, we can continue as below:

where $l_{1} = \left(18a^{2}M_{1}^{4}+8b^{2}M_{1}^{2}\right).$

If we have $u, u_{1}\in R^{2}$ and $t\in\lbrack0, T],$

Finally, if we have $w, w_{1}\in R^{2}$ and $t\in\lbrack0, T],$

where $l_{3} = \frac{3}{2}r^{2}.$

So, Condition (ⅰ) is satisfied easily.

Now we will verify the second condition (ⅱ) for our system.

$\forall(t, v(t))\in R^{2}\times\lbrack0, T]$; then, we will show that

where

and while satisfying the condition

Now, we continue with the second equation.

$\forall(t, u(t))\in R^{2}\times\lbrack0, T]$; then, we will show that

where

and while satisfying the condition

Finally, let us take the last equation, $\forall(t, w(t))\in R^{2}\times \lbrack0, T]$; then, we will show that

where

and satisfying the condition

So, ıf the conditions below are satisfied, the model has a unique solution.

Under the above conditions, we can conclude that our system admits a unique exact system of positive solutions.

3.1.   Existence and uniqueness for integral expressions

We now consider the model as below;

We now apply the integration to both sides to get

Let us consider the following mapping

We evaluate

So we have

We next evaluate

Similarly,

We now consider the model with the Caputo derivative.

If we convert the system above into a fractional inrtegral equation, we get

Again,

To proceed, we use the Hölder inequality to obtain

Following the procedure presented above, we concluded that

We now evaluate the Lipschitz condition for all submappings:

We make use of the Hölder inequality to establish the Lipschitz

Thus,

with the same procedure, we have

The existence and uniqueness of the model with the Caputo-Fabrizio and the Atangana-Baleanu can be obtained similarly. In the next section, we present a numerical solution and apply it to solve the model presented here.

3.2.   Numerical scheme for model

The existence and uniqueness of the model with the Caputo-Fabrizio and the Atangana-Baleanu can be obtained similarly. In the next section, we present a numerical solution and apply it to solve the model presented here.

Case 1: Let us consider the following general Cauchy problem:

where $f\left(.\right)$ is twice differentiable and verifies the conditions presented above into the integral equation as follows:

We evaluate the integral at $t = t_{n+1}$:

within $[t_{2}, t_{n+1}]$ we approximate $f(\tau, y(\tau))$ with Newton's interpolation:

Replacing all into the main equation, we obtain

Thus,

We must check the stability of the approximation. To achieve this, we consider the perturbation terms $\widetilde{y}_{0}, \widetilde{y}_{1}$ and $\widetilde{y} _{j}$ for $j = 2, ..., N.$

Taking the difference between the perturbed approximation and the approximation yields

Using the Lipschitz condition of $f(., y)$ yields

Let $\overline{\alpha} = \underset{0\leq j\leq N}{\max}\left\{ \left\vert \widetilde{y}_{j}\right\vert \right\};$ then,

Case 2: We consider the case where the Cauchy problem is with the Caputo derivative.

From the above, we can build the following:

at $t = t_{n+1}$

We have that, within $[0, t_{1})$ we can approximate $f(\tau, y(\tau))$ as

and in $[t_{1}, t_{2}]$ we approximate $f(\tau, y(\tau))$ as

Then, replacing in their respective integral yields

Thus,

We evaluate first

Therefore,

Therefore, we have

Case 3: We consider the case where the Cauchy problem has the Atangana-Baleanu derivative

Applying the Atangana-Baleanu fractional integral yields

At $t = t_{n+1},$ we have

Then, we have

where

3.3.   Numerical simulations

As indicated earlier the biological role can be described by the last equation to capture different dynamical patterns found in a real-world situation; we should also note that the classical differential operators used to model this process is based on the Delta-Dirac kernel. Of course, this is because the first derivative is the convolution of itself and the Delta-Dirac function. This gives such operators fewer properties to capture more patterns or nonlocal patterns that are found in real-world problems. On the other hand differential operators based on the power-law kernel have been recognized as a good mathematical tool to replicate processes with power-law behaviors. Thus, to include power-law processes into the existing mathematical model, the classical differential operator will be replaced by the Caputo-power law derivative and the derivative based on the generalized Mittag-Leffler kernel as this kernel provides a crossover behavior from the stretched exponential to the power law, which is a relaxation that is found in many biological processes.

where the initial condition is considered as

In Figures 1–3, we show the numerical simulations with the parameters $a = 0.1$, $b = 3,$ $c = 1,$ $d = 5$, $I = 6,$ $s = 40$ $r = 0.001,$ $v_{rest} = -1.6$ and $\alpha = 1.$

In Figures 4–6, we show the numerical simulations with the parameters $a = 1$, $b = 3,$ $c = 1,$ $d = 5$, $I = 6,$ $s = 40$ $r = 0.001,$ $v_{rest} = -1.6$ and $\alpha = 1.$

In Figures 7–9, we show the numerical simulations with the parameters $a = 2$, $b = 3,$ $c = 1,$ $d = 5$, $I = 2,$ $s = 4$ $r = 0.0001,$ $v_{rest} = -1.6$ and $\alpha = 1.$

In Figures 10–12, we show the numerical simulations with the parameters $a = 2$, $b = 3,$ $c = 1,$ $d = 5$, $I = 2,$ $s = 4$ $r = 0.0001,$ $v_{rest} = 10$ and $\alpha = 1.$

Comments for figures: As previously mentioned, nonlocal operators are based on the convolution of the rate of change and a few significant mathematical functions that share characteristics with observations from the real world. The figures that were produced by the model using the generalized Mittag-Leffler function show this. These graphs show the relationship between membrane potential, the adaption current and measurements of the transit of sodium and potassium ions through fast ion channels. We have a parametric representation of these functions in Figures 1 through 12. The Atangana-Baleanu fractional derivative, which is based on the generalized Mittag-Leffler function, intensifies the chaotic behaviors of the spiking-bursting conduct of the membrane potential in these pictures. In particular, the generalized Mittag-Leffler function-induced trend with crossover from the stretched exponential to the power law may be seen. This is distinct from the results that were achieved when the model was built using a classical differential operator.

4.   Conclusions

We considered the Hindmarsh-Rose model to further investigate the dynamic process of the spiking-bursting conduct of the membrane potential seen in experiments finished with a single neuron. The model is composed of a set of three nonlinear ordinary differential equations, where the function $v(t)$ represents the membrane potential and the function $u(t)$ represents the measurement of the movement of sodium and potassium ions through fast ion channels. The adaptation current is represented by the function $w(t)$. There were some theoretical evaluations offered. We took into account a few common Cauchy issues and put a modified plan into practice using Newton's polynomial interpolation. In order to show how nonlocal behaviors will affect the model, we presented its consistency and convergence analyses.