




































©2023 Ada Academica https://adac.eeEur. J. Math. Anal. 3 (2023) 2doi: 10.28924/ada/ma.3.2
Analysis of Neuronal Oscillations of Fractional-Order Morris-Lecar Model

Tahmineh Azizi
Department of Mechanical Engineering, Florida State University, USA

Correspondence: tazizi@fsu.edu

Abstract. Fractional calculus is a new approach for modeling biological and physical phenomena withmemory effects. Fractional calculus uses differential and integral operators including non-integerorders to study the non-linear behavior of physical and biological systems with some degrees offractionality or fractality. Since the long memory properties of neuronal responses can be betterexplained using fractional derivative, in this study we generalize the integer-order Morris-Lecarmodel in the fractional-order domain to better modeling of neuron dynamics. To investigate thecomplex spiking patterns of fractional-order Morris-Lecar neural system the fractional calculus hasbeen applied to build this new mathematical model. We compare the results with integer-order Morris-Lecar model. The analytical solutions of these equations cannot explicitly be obtained. Therefore, tofind the dynamical behaviors of solutions, we used approximation and numerical schemes. Dependingon the different parameters values for 0 < η ≤ 1, the fractional-order Morris-Lecar reproducesquiescent, spiking and bursting activities the same as its original model but for higher input current.We numerically discover the hopf bifurcation, saddle node bifurcation of limit cycle and homoclinicbifurcation for this model for different input current and derivative orders. Taking the advantages ofthe fractional order derivative, for a variety of orders, we define different classes of this model whichhelps to better extract all the complicated dynamics of this single neuron model.

1. Introduction
Recently, fractional calculus has been frequently used by many researchers in biology, physics,chemistry and biochemistry, hydrology, medicine, and finance and its application in modelingcomplex phenomena has increased its popularity and the number of publications in above area [1–6].The main characteristic of fractional order derivative is called the "memory effect" and it has beenexperimentally proved that the fractional order differential or integral equations models are morerealistic to demonstrate the complex behavior of some biological or physical systems includingfractality and memory compared to their ODEs of integer-order systems [2, 6]. Complexity in thiscontext combined the recent advances in neuroscience with the concepts from fractal geometry andnonlinear dynamics to form a new approach within the life sciences which is useful in controlling
Received: 17 Apr 2022.
Key words and phrases. Fractional calculus; NSFD method; Hopf bifurcation; Homoclinic bifurcation; Grunwald-Letinkov method. 1

https://adac.ee
https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 2the dynamics of fractal processes in this area [7, 8].Understanding the complicated functioning of the neuronal cells and exploring the molecular andcellular mechanisms of their network have been one of the greatest challenges in different fieldsof science . The progress and advances in computational neuroscience could help scientists tobetter understanding of the performance of brain and neuron cells and better fighting with diseasesrelated to neuron cells such as Parkinson’s and depression.Non-linear dynamical system theory has a very important role in the computational neuroscienceresearch [9–13]. In 1948 Hodgkin by injecting a dc-current of varying amplitude discovered thatsome preparations could show repetitive spiking activities with arbitrarily low frequencies, whilethe others discharged in a narrow frequency band [9, 13–15]. His finding motivated Rinzel andErmentrout to discover that different bifurcation mechanisms of excitability may cause the differencein neuronal behavior [9,16,17]. Basically, if assume the applied current Iapp as a control parameter,we can easily see the transition in behavior of a neuron which corresponds to a bifurcation fromequilibrium to a limit cycle attractor. That is when Iapp is small, the cell remains quiescent andwith increasing the injected current, the cell starts to fire repetitive spikes [9–13,18,19].According to Moaddy. K, et. al [20], the fractional-order models can better explain the long memorydependence of the neuron response. One of the most interesting properties of neural system isadaptation to changes in stimulus. It has been shown that a single neuron has a single time scaleadaptation, however, there are some neurons with multiple time scale adaptation to responsesconsistent with fractional-order derivatives, means that the firing rate for these neurons acts asfractional derivative of slowly varying stimulus parameters [21,22]. Therefore, the other advantagesof neuronal fractional derivatives is their ability to adapt to changes in stimulus in different timescales. Moreover, Shi. M, et. al [23] proved that the fractional-order derivative demonstrates thereal dielectric behaviors and the history memory property of membranes, cells and so on. Accordingto their finding, non integer derivative activates the slow ion channel with higher speed, and helpsto activate fast spiking modulation which forms different kinds of bursting behaviors. One importantfact about using fractional order model their ability to display different dynamical behaviors such aschaotic and periodic firing for the same parameter values as the order of derivative is varying [24].On the other hand, using these new parameters as the order of fractional derivative operatorsenhances the controllability of behavior of the neuron cells [25]. The fractional order models playan important role in determining the firing properties of neuronal models and these importantproperties of models with fractional order derivatives including depiction of long term memory andthe multiple time scale adaptation have motivated many efforts to use them as a perfect frameworkto cover the complicated dynamics of many different neuronal cells such as fractional cable model,Izhikevich neuron model, and FitzHugh-Rinzel bursting neuron model [26–28].Due to the complexity of nerve systems, it is impossible to fully understand the various phenomena

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 3in neuroscience only using the integer order models, since it does not meet all neuronal propertiesand complicated behaviors of neurons. Thus, in this study we use the fractional calculus to explorethe dynamics of fractal processes using the fractional calculus and apply this dynamical approach onMorris-Lecar model to catch all the spiking properties of this neuron and to simulate the fluctuationsof this single neuron cell and obtain biological physiological characteristics of it. We have selectedMorris-Lecar model because it is a reduced and simpler version of the Hodgkin-Huxley equationsand preserves many important characteristics of neuronal dynamics such as generation of actionpotentials, threshold for firing spike, and sustained oscillations with increasing the applied current.Because the solutions of fractional Morris Lecar model (FML) may not be explicitly obtained, weuse numerical methods to approximate the solutions of this model. We compare these results withits original integer order model using phase portrait analysis. By considering this fractional ordermodel, we can explain all the possible geometric mechanisms underlying each of neuronal activitiesof Morris-Lecar model.
2. Grünwald-Letinkov approximation

We define the fractional differential as the following form [29,30]
DγY (t) = f (t, Y (t)), Y (t0) = Y0

where γ > 0 represents the order of derivative and Dγ denotes the fractional derivative which isgiven by:
DγY (t) = Jk−γDkY (t)

where γ ∈ (k − 1, k ], for k = 1, 2, . . . and integral operator Jk called the Riemann-Liouville of
kth-order which is obtained by the following formula

JkY (t) =
1

Γ(k)

∫ t

0

(t − τ)(k−1)Y (τ) dτ, t > 0

where Γ(.) denotes the gamma function.To apply the Micken’s (NSFD) [31–33], we need to find the fractional order derivative using theGrünwald- Letinkov (G-L) approximation for model equations as the form
DγY (t) = lim

s→0
s−γ

T∑
i=0

(−1)i
(
γ

i

)
Y (t − i s) (1)

where T = [t]/s and [.] used to show the integer value and s represents the step size. Thus,equation (1) would be discretized as
T∑
i=0

Cγi Y (tk−i) = f (tk , Y (tk)) (2)

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 4where tk = k s and Cγi are the coefficients for (G-L) approximation written as
Cγi =

[
i − 1− γ

i

]
Cγi−1, Cγ0 = s−γ i = 1, 2, . . .

Next we introduce the non-standard finite difference schemes.To discretize a systems of differential equations, both ordinary differential equations (ODEs) andpartial differential equations (PDEs), one may apply the Mickens NSFD discretization methodwhich is more flexible in construction rather than standard finite difference method and thereforehas better performance. This method checks the positivity of solutions and is concerned aboutboundedness and monotonicity of them. Another advantage of using NSFD schemes is their abilityto preserve the structure and properties of the systems of differential equations and therefore, weapply NSFD schemes on the general compartmental model in the form:
d Y

dt
= f (Y ) (3)

However, to use the non-standard scheme we need to check that if non-local approximation is usedand or we need to have a non traditional discretization of derivatives and also we may need to use anon-negative function Φ(h) = s+O(s2). To apply NSFD scheme, we consider a grid tk = t0+k s ,such that s > 0, and we approximately write the discretized function Y as Yk ≈ Y (tk). Next, wediscretize (3):
d Y

dt
=
Yk+1 − Yk

Φ(s)
+O(Φ(s)) (4)

when s → 0 we have
d Y

dt
≈
Yk+1 − Yk

Φ(s)
(5)

where real valued Φ(s) as a function of the step size s need to satisfy the following properties [34]:(I) Φ(s) = s +O(s2),(II) Φ(s) ∈ (0, 1), ∀ s ∈ (0,∞)Here, the equality (4) is equivalent with the integer order derivative as follow:
d Y

dt
= lim
s→0

[
Y (t + s)− Y (t)

Φ(s)
+O(Φ(s))

]
= lim
s→0

[
Y (t + s)− Y (t)

s

]
lim
s→0

[
s

Φ(s)

]
+ lim
s→0

O(Φ(s))

= Ẏ (t)

As s → 0 the discrete form in (4) converges to its associated continuous derivative. NSFD methodsare convergent without any restriction related to step size s but this is not always true for SFDmethods which depend on the step size s . Moreover, when we discretize a system using NSFD

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 5method, if the original system is persistent, and solutions are stable and convergent, these proper-ties remain the same after discretization, but not for the case we use SFD to discretize the systemof differential equations.
3. Description of model equations

To demonstrate the generation of action potential, Kathleen Morris and Harold Lecar proposeda simple model, Morris-Lecar model, in 1981 [35] that is a reduction version of the four dimensionalHodgkin-Huxley model preserving the main properties of spike generations with much simplermathematical and computational analysis [35, 36]. This model describes the electrical activities ofneurons using a system of non-linear ordinary differential equations and includes three channelsa potassium channel, a leak and a calcium channel and has the following form
CM

dV

dt
= Iapp − gL(V − EL)− gKw(V − EK)− gCam∞(V )(V − ECa) = Iapp − Iion(V, w),

dw

dt
= φ(n∞(V )− w)/τw (V ),

(6)
where

m∞(V ) =
1

2
[1 + tanh((V − V1)/V2)], (7)

τn(V ) = 1/cosh((V − V3)/(2V4)), (8)
n∞(V ) =

1

2
[1 + tanh((V − V3)/V4)]. (9)

and
Iion(V, w) = gL(V − EL) + gKw(V − EK) + gCam∞(V )(V − ECa) (10)

where V demonstrates membrane potential, and w the activation variable of the persistent K+current, so it is a two-dimensional vector (V, w). EK , ECa, and EL denote the Nernst equilibriumpotentials. Iapp demonstrates the injected current and Iion the ionic current. Parameter φ is atemperature factor. gL is leak membrane conductance, gK is potassium membrane conductance and
gCa is calcium membrane conductance. Moreover, CM is the total membrane capacitance. Also, thevoltage-sensitive steady-state activation function m∞(V ) and n∞(V ), and the time constant τw (V )can be measured experimentally. The non-linear dynamics of the original Morris-Lecar model havebeen studied by different researchers during recent decades [18,37–43]. In the next section we willlook at the fractional-order Morris-Lecar model and its spiking patterns.
3.1. Fractional Morris-Lecar model. Now, we apply the basic theorems of the fractional calculuson model (6). In Morris-Lecar model, we write the total membrane current to being the sum of ioniccurrents and the capacitive current:

Iapp = Iion + ICM

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 6For the fractional-order model, we define:
ICM = CM Dη V

where, 0 < η ≤ 1 and Dη is defined in the following form [44]:
Dη V (t) = lim

h→0
h−η

[t]
h∑
i=0

(−1)i
(
η

i

)
V (t − i h) (11)

We do the same for the second equation:
Dη w(t) = lim

h→0
h−η

[t]
h∑
i=0

(−1)i
(
η

i

)
w(t − i h) (12)

where [t] denotes the integer part of t and h is the step size.After discretization, (11) and (12) become:
∑ [t]

h

i=0 C
η
i V (tk−i) = f (tk , V (tk)),

∑ [t]
h

i=0 C
η
i w(tk−i) = g(tk , w(tk)),

where tk = kh for k = 1, 2, 3, . . . and Cηi are the Grunwald-Letinkov coefficients as:
Cηi =

(
1−

1 + η

i

)
Cηi−1, Cη0 = h−η i = 1, 2, . . .

Then, we apply the non-standard finite difference (NSFD) schemes proposed by Mickens [31–33]and replace the step size h by a function ψ(h). Next, we discretize the equations (11) and (12)following the Grunwald-Letinkov discretization, using V (tk) = Vk , w(tk) = wk , we have:
CM

∑k+1
i=0 C

η
i Vk+1−i = Iapp − gL(Vk − EL)− gKwk(Vk − EK)− gCam∞(Vk)(Vk − ECa),

∑k+1
i=0 C

η
i wk+1−i = φ(n∞(Vk)− wk)/τwk (Vk),

(13)
where

m∞(Vk) =
1

2
[1 + tanh((Vk − V1)/V2)], (14)

τn(Vk) = 1/cosh((Vk − V3)/(2V4)), (15)
n∞(Vk) =

1

2
[1 + tanh((Vk − V3)/V4)]. (16)

and
Iion(Vk , wk) = gL(Vk − EL) + gKwk(Vk − EK) + gCam∞(Vk)(Vk − ECa) (17)

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 7

Table 1. Parameter values for the fractional-order Morris-Lecar model.
After some algebra,

Vk+1 =
Iapp −

∑k+1
i=1 C

η
i Vk+1−i − gL(Vk − EL)− gKwk(Vk − EK)− gCam∞(Vk)(Vk − ECa)

CM C
η
0

wk+1 =
φ(n∞(Vk)− wk)− τwk (Vk)

∑k+1
i=1 C

η
i wk+1−i

Cη0 τwk (Vk) (18)
where

Cη0 = ψ(h)−η, ψ(h) = sin(h)

The fractional-order Morris-Lecar model displays different ranges of dynamics such as hopfbifurcation, saddle node on invariant limit cycles (SNLC) and homoclinic bifurcation. We keepthe same biological parameters as the original Morris-Lecar model. We have represented theseparameters for these three different dynamics in table (1) [18]. We assume Iapp as a controlparameter for numerical simulations.
3.2. Local stability analysis of fractional order Morris-Lecar model. In neuroscience, it’s usuallyhard to extract analytically the dynamics of the neuronal systems and we may need to use somegeometrical and qualitative techniques such as phase portrait analysis. Phase portraits demonstratethe evolution of state variables in time with different initial states. By looking at the phase portrait,we can observe the qualitative behavior of the system without knowing the model equations. Toanalyze the local dynamics of fractional order Morris-Lecar model, we apply a useful theoremin dynamical systems theory, called the Hartman-Grobman theorem [45–47]. According to this

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 8theorem non-linear fractional order Morris-Lecar system{
Vk+1 = F (Vk , wk)

wk+1 = G(Vk , wk)
(19)

sufficiently near equilibrium (V, w) = (V ∗, w∗) is locally topologically equivalent to the linear partof the system. First we transform the fixed point (V ∗, w∗) of the system (19) to the origin by thetranslation V = V ∗ + V̄ and w = w∗ + w̄ . If we split off the linear part of the system from itsnon-linear part, we have
[
V̄

w̄

]
=


∂F (V̄ , w̄)

∂V̄

∂F (V̄ , w̄)

∂w̄

∂G(V̄ , w̄)

∂V̄

∂G(V̄ , w̄)

∂w̄


[
V̄

w̄

]
+

F̃ (V̄ , w̄)

G̃(V̄ , w̄)

 (20)
where F̃ and G̃ represent the non-linear part of the system (19) and

∂F (V̄ , w̄)

∂V̄
=
−
∑k+1
i=1 C

η
i − gL − gKw̄ − gCam∞(V̄ )

CM C
η
0

≡ A

∂F (V̄ , w̄)

∂w̄
=
−gK V̄
CM C

η
0

≡ B

∂G(V̄ , w̄)

∂V̄
=

φn′∞(V̄ )

Cη0 τwk (V̄ )
≡ C

∂G(V̄ , w̄)

∂w̄
=
−φ− τwk (V̄ )

∑k+1
i=1 C

η
i

Cη0 τwk (V̄ )
≡ D

(21)

To find the stability of the interior equilibrium point of the system, we need to look at the linearpart of (20) which is given by (21). At first, we assume that for the voltage Ek < V̄ < ECa. Thenwe have B < 0, C > 0, and D < 0 and A can be either positive or negative. m∞(V̄ ) which definesthe slope of the calcium activation function can make A > 0. On the other hand, for A < 0, theequilibrium point is asymptotically stable because A + D < 0 and AD − BC > 0. Moreover, for
B < 0 the equilibrium point is stable because the negativity of slope of V̄ -nullcline, −AB < 0.For the case that A > 0 the equilibrium point is a saddle point and unstable because of thepositivity of slope of the w̄-nullcline −CD > 0.For −AB > −C

D , the equilibrium point is a saddle point and unstable because AD − BC < 0.However, for −AB < −C
D and A + D < 0, the equilibrium point is stable and for −AB < −C

D and
A+D > 0, the equilibrium point is unstable.Finally, for the case that A > 0 and the speed of potassium dynamics φ is small, then the equilibriumpoint is unstable.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 94. Numerical results
In this section, we use some numerical simulations to study the qualitative behavior such aslocal bifurcations of the fractional order Morris-Lecar model for different fractional order η andapplied current Iapp . One of the most common types of bifurcation in neuroscience is saddle nodebifurcation of limit cycle or SNLC, and this bifurcation occurs when with increasing the controlparameter, here, applied current Iapp , two stable and unstable limit cycles which are associated tothe stable node and saddle point respectively, close to each other, collide and at the bifurcationtime, a limit cycle appears. With increasing Iapp further, this limit cycle disappears. Figures(1)-(4) exhibit different spiking behaviors for fractional Morris-Lecar model (13), when we increase

Iapp = 5, 30, 45, 100 using SNLC parameters value in table (1). The solution of the integer ordermodel (6) has been demonstrated in the third row to compare with fractional-order Morris-Lecarmodel of different order.
For the case of hopf bifurcation in figures (5)-(9), with increasing the applied current Iapp , themodel (13) displays the occurrence of limit cycle corresponding to hopf bifurcation like the originalmodel (6) but for orders η = 0.3, .0.5, 0.7, 0.9 the fractional order model needs greater value forinput current to start the bifurcation.The topological normal form of the model (13) in polar coordinate for the case of hopf bifurcationhas the form:


∑k+1
i=0 C

η
i rk+1−i = αr + a r3

∑k+1
i=0 C

η
i θk+1−i = ω0 + βr2

(22)
After simplification, the fractional-order system which is linear and time-invariant has the followingform: 

rk+1 =
αr + a r3 −

∑k+1
i=1 C

η
i rk+1−i

Cη0

θk+1 =
ω0 + βr2 −

∑k+1
i=1 C

η
i θk+1−i

Cη0

(23)
where, α and ω0 represent the real part and imaginary part of the eigenvalues of the jacobian matrixfor the model (13) around its equilibrium point respectively, a is called first Lyapunov coefficient.For a > 0 there should exist an unstable limit cycle, bifurcating from the equilibrium and it indicatesthe appearance of subcritical hopf bifurcation and for a < 0 we have stable limit cycle solutionand supercritical hopf bifurcates from the equilibrium. Here, β does not have any dynamical effect.
θ represents the angle of oscillations. If θ̇ > 0, it means the frequency of damped or sustained

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 10

Figure 1. Occurrence of saddle node bifurcation of limit cycle or SNLC in FractionalMorris-Lecar model (13) for Iapp = 5, third row displays the trajectory of the originalmodel (6) with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 11

Figure 2. Occurrence of saddle node bifurcation of limit cycle or SNLC in FractionalMorris-Lecar model (13) for Iapp = 30, third row displays the trajectory of theoriginal model (6) with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 12

Figure 3. Occurrence of saddle node bifurcation of limit cycle or SNLC in FractionalMorris-Lecar model (13) for Iapp = 45, third row displays the trajectory of theoriginal model (6) with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 13

Figure 4. Occurrence of saddle node bifurcation of limit cycle or SNLC in FractionalMorris-Lecar model (13) for Iapp = 100, third row displays the trajectory of theoriginal model (6) with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 14

Figure 5. Occurrence of hopf bifurcation in Fractional Morris-Lecar model (13) for
Iapp = 20, third row displays the trajectory of the original model (6) with the sameapplied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 15

Figure 6. Occurrence of hopf bifurcation in Fractional Morris-Lecar model (13) for
Iapp = 88, third row displays the trajectory of the original model (6) with the sameapplied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 16

Figure 7. Occurrence of hopf bifurcation in Fractional Morris-Lecar model (13) for
Iapp = 90, third row displays the trajectory of the original model (6) with the sameapplied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 17

Figure 8. Occurrence of hopf bifurcation in Fractional Morris-Lecar model (13) for
Iapp = 95, third row displays the trajectory of the original model (6) with the sameapplied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 18

Figure 9. Occurrence of hopf bifurcation in Fractional Morris-Lecar model (13) for
Iapp = 220, third row displays the trajectory of the original model (6) with the sameapplied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 19oscillations around ω0 is increasing. On the other hand, for θ̇ < 0 the frequency of damped orsustained oscillations around ω0 is decreasing.In neuroscience point of view, the hopf bifurcation happens when the behaviors of neuron changefrom resting to spiking (the stable constant solutions are corresponding to the resting state andspiking state shows the existence of periodic solutions).The other common type of dynamical behavior for a neuron cell occurs when with increasing thecontrol parameter, a saddle point and a limit cycle collide, this bifurcation called saddle-homoclinicbifurcation. The period of the periodic orbit that appears at the moment of bifurcation goes to infinityand with further increasing of control parameter this periodic orbit disappears. Figures (10)-(14),demonstrate the appearance and disappearance of saddle-homoclinic bifurcation in the model (13)with increasing the applied current Iapp = 23, 40, 50, 60, 70 like the original model (6) but liketwo previous bifurcations, for fractional order model of orders η = 0.3, .0.5, 0.7, 0.9 the neuronneeds higher input current Iapp to bifurcate.In neuroscience point of view, when saddle homoclinic bifurcation happens, we expect theappearance or disappearance of spiking behavior.
5. Discussion

Fractional-order excitable systems can be physically considered as a memory dependent phe-nomenon which display oscillatory behaviors for certain types of neuron models. In this research,we have studied the neuronal spiking patterns of fractional Morris-Lecal neuron model where thefractional-orders could change the responses of the model from periodic to non-periodic, and wehave compared its dynamics to the original Morris-Lecar model. The original Morris-Lecal neuronmodel which is a reduction version of Hodgkin-Huxley model includes two equations with integerorder derivatives and three ionic channels, a potassium channel, a leak and a calcium channel.We have preserved the same ionic channels and to find the fractional order Morris-Lecar model,we have applied the non-standard finite difference (NSFD) schemes on this system of equationssince they have a better performance than standard finite difference methods. Then we have dis-cretized the model using the Grunwald-Letinkov discretization. We use effective numerical methodsto display the solution of fractional order Morris-Lecar model. To explore the exciting behaviorsof fractional Morris-Lecar model, we have conducted different numerical simulations with changinginput currents. It was obvious that the solutions depends on the fractional-order parameters. Wehave shown that the fractional Morris-Lecar model with the same biological parameters values asthe original model, displays the same firing patterns such as quiescent and spiking behaviors butfor different values of input currents. We have noticed that in this case, the saddle node bifurcationof limit cycle (SNLC), hopf bifurcation and saddle-homoclinic bifurcation happen at larger values forinjected current compare to the original model and we have derived these bifurcations analyticallyusing rigorous normal form theory. Similar to the original Morris-Lecar model, its fractional-order

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 20

Figure 10. Occurrence of saddle-homoclinic bifurcation in Fractional Morris-Lecarmodel (13) for Iapp = 23, third row displays the trajectory of the original model (6)with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 21

Figure 11. Occurrence of saddle-homoclinic bifurcation in Fractional Morris-Lecarmodel (13) for Iapp = 40, third row displays the trajectory of the original model (6)with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 22

Figure 12. Occurrence of saddle-homoclinic bifurcation in Fractional Morris-Lecarmodel (13) for Iapp = 50, third row displays the trajectory of the original model (6)with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 23

Figure 13. Occurrence of saddle-homoclinic bifurcation in Fractional Morris-Lecarmodel (13) for Iapp = 60, third row displays the trajectory of the original model (6)with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 24

Figure 14. Occurrence of saddle-homoclinic bifurcation in Fractional Morris-Lecarmodel (13) for Iapp = 70, third row displays the trajectory of the original model (6)with the same applied current.

https://doi.org/10.28924/ada/ma.3.2


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 25model undergoes a transition between integrator and resonator. When saddle-node bifurcationhappens, the neuron is called an integrator means that there is no damped subthreshold oscilla-tions. On the other hand, when Hopf bifurcation happens the neuron is called a resonator withdamped subthreshold oscillations. Using the fractional order derivative, we have added a newparameter as the order of derivatives that helped us to control the spiking patterns of the neuroncell. Taking the advantages of this type modeling, we investigated how the classical order systemschanges its complex dynamics such as firing patterns and also firing frequency, when they turn tobe fractional order systems. This work improved the preceding ones [37,38] by discovering differentattractors of the system for different fractional orders and keeping the same biological parameters.As a result, fractional order plays a key role in describing the firing patterns and characterizing thememory effect of neurons which helps to control the long term dependency of the neuron responsesby adding extra freedom to the system. Moreover, using this fractional operator could display morenon-local natural dynamics as a sign of fractal behaviors compared to the integer order model. Thisdifferentiation operator which is a combination of fractal and fractional differentiation indicated theimportance of using fractal geometry to study the neural dynamic systems.
References

[1] R. Hilfer, Applications of fractional calculus in physics, World Scientific, Singapore, 2000. https://doi.org/10.
1142/3779.[2] F.A. Rihan, D. Baleanu, S. Lakshmanan, R. Rakkiyappan, On fractional SIRC model with salmonella bacterial infection,2014 (2014) 136263. https://doi.org/10.1155/2014/136263.[3] F.A. Rihan, S. Lakshmanan, A.H. Hashish, R. Rakkiyappan, E. Ahmed, Fractional-order delayed predator–prey sys-tems with Holling type-II functional response, Nonlinear Dyn. 80 (2015) 777–789. https://doi.org/10.1007/
s11071-015-1905-8.[4] F.A. Rihan, A. Hashish, F. Al-Maskari, et al. Dynamics of tumor-immune system with fractional-order, J. Tumor Res.2 (2016) 109.[5] J. Gómez-Aguilar, M. López-López, V. Alvarado-Martínez, D. Baleanu, H. Khan, Chaos in a cancer model via fractionalderivatives with exponential decay and Mittag-Leffler law, Entropy. 19 (2017) 681. https://doi.org/10.3390/
e19120681.[6] M. Zeinadini, M. Namjoo, Approximation of fractional-order Chemostat model with nonstandard finite differencescheme, Hacettepe J. Math. Stat. 46 (2017) 469-482.[7] D. Matignon, Stability results for fractional differential equations with applications to control processing, Comput.Eng. Syst. Appl. 2 (1996) 963-968.[8] E. Ahmed, A.M.A. El-Sayed, H.A.A. El-Saka, On some Routh–Hurwitz conditions for fractional order differentialequations and their applications in Lorenz, Rössler, Chua and Chen systems, Phys. Lett. A. 358 (2006) 1–4. https:
//doi.org/10.1016/j.physleta.2006.04.087.[9] E.M. Izhikevich, Dynamical systems in neuroscience, MIT Press, 2007.[10] E.M. Izhikevich, Simple model of spiking neurons, IEEE Trans. Neural Netw. 14 (2003) 1569–1572. https://doi.
org/10.1109/tnn.2003.820440.[11] E.M. Izhikevich, Which Model to Use for Cortical Spiking Neurons?, IEEE Trans. Neural Netw. 15 (2004) 1063–1070.
https://doi.org/10.1109/tnn.2004.832719.

https://doi.org/10.28924/ada/ma.3.2
https://doi.org/10.1142/3779
https://doi.org/10.1142/3779
https://doi.org/10.1155/2014/136263
https://doi.org/10.1007/s11071-015-1905-8
https://doi.org/10.1007/s11071-015-1905-8
https://doi.org/10.3390/e19120681
https://doi.org/10.3390/e19120681
https://doi.org/10.1016/j.physleta.2006.04.087
https://doi.org/10.1016/j.physleta.2006.04.087
https://doi.org/10.1109/tnn.2003.820440
https://doi.org/10.1109/tnn.2003.820440
https://doi.org/10.1109/tnn.2004.832719


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 26

[12] E.M. Izhikevich, Neural excitability, spiking and bursting, Int. J. Bifurcation Chaos. 10 (2000) 1171–1266. https:
//doi.org/10.1142/s0218127400000840.[13] F.C. Hoppensteadt, E.M. Izhikevich, Weakly Connected Neural Networks, Springer New York, 2012. https://doi.
org/10.1007/978-1-4612-1828-9.[14] A.L. Hodgkin, The local electric changes associated with repetitive action in a non-medullated axon, J. Physiol. 107(1948) 165–181. https://doi.org/10.1113/jphysiol.1948.sp004260.[15] A.L. Hodgkin, A.F. Huxley, A quantitative description of membrane current and its application to conduction andexcitation in nerve, J. Physiol. 117 (1952), 500–544. https://doi.org/10.1113/jphysiol.1952.sp004764.[16] J. Rinzel, G.B. Ermentrout, Analysis of neural excitability and oscillations, MIT Press, Cambridge, MA, 1989.[17] J. Rinzel, G.B. Ermentrout, Analysis of neural excitability and oscillations, Methods in neuronal modeling, MITpress Cambridge, MA, 1998.[18] G.B. Ermentrout, D.H. Terman, Mathematical Foundations of Neuroscience, Springer New York, 2010. https:
//doi.org/10.1007/978-0-387-87708-2.[19] B. Ermentrout, Type I membranes, phase resetting curves, and synchrony, Neural Comput. 8 (1996) 979–1001.
https://doi.org/10.1162/neco.1996.8.5.979.[20] K. Moaddy, A.G. Radwan, K.N. Salama, S. Momani, I. Hashim, The fractional-order modeling and synchronizationof electrically coupled neuron systems, Computers Math. Appl. 64 (2012) 3329–3339. https://doi.org/10.1016/
j.camwa.2012.01.005.[21] B.N. Lundstrom, M.H. Higgs, W.J. Spain, A.L. Fairhall, Fractional differentiation by neocortical pyramidal neurons,Nat. Neurosci. 11 (2008) 1335–1342. https://doi.org/10.1038/nn.2212.[22] S.A. Malik, A.H. Mir, FPGA realization of fractional order neuron, Appl. Math. Model. 81 (2020) 372–385. https:
//doi.org/10.1016/j.apm.2019.12.008.[23] M. Shi, Z. Wang, Abundant bursting patterns of a fractional-order Morris–Lecar neuron model, Commun. NonlinearSci. Numer. Simul. 19 (2014) 1956–1969. https://doi.org/10.1016/j.cnsns.2013.10.032.[24] D. Jun, Z. Guang-jun, X. Yong, Y. Hong, W. Jue, Dynamic behavior analysis of fractional-order Hindmarsh–Roseneuronal model, Cogn. Neurodyn. 8 (2013) 167–175. https://doi.org/10.1007/s11571-013-9273-x.[25] Z. Wang, X. Wang, Y. Li, X. Huang, Stability and hopf bifurcation of fractional-order complex-valued sin-gle neuron model with time delay, Int. J. Bifurcation Chaos. 27 (2017) 1750209. https://doi.org/10.1142/
s0218127417502091.[26] B.I. Henry, T.A.M. Langlands, S.L. Wearne, Fractional Cable Models for Spiny Neuronal Dendrites, Phys. Rev. Lett.100 (2008) 128103. https://doi.org/10.1103/physrevlett.100.128103.[27] M.F. Tolba, A.H. Elsafty, M. Armanyos, L.A. Said, A.H. Madian, A.G. Radwan, Synchronization and FPGA realizationof fractional-order Izhikevich neuron model, Microelectronics J. 89 (2019) 56–69. https://doi.org/10.1016/j.
mejo.2019.05.003.[28] A. Mondal, S.K. Sharma, R.K. Upadhyay, A. Mondal, Firing activities of a fractional-order FitzHugh-Rinzelbursting neuron model and its coupled dynamics, Sci. Rep. 9 (2019) 15721. https://doi.org/10.1038/

s41598-019-52061-4.[29] K.S. Miller, B. Ross, An introduction to the fractional calculus and fractional differential equations, Wiley, 1993.[30] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equa-tions, to methods of their solution and some of their applications, Elsevier, 1998.[31] R.E. Mickens, Nonstandard finite difference models of differential equations, World Scientific, 1994.[32] R.E. Mickens, Nonstandard finite difference schemes for reaction-diffusion equations, Numer. Meth. Part. Differ. Equ.15 (1999) 201–214. https://doi.org/10.1002/(sici)1098-2426(199903)15:2<201::aid-num5>3.0.co;2-h.

https://doi.org/10.28924/ada/ma.3.2
https://doi.org/10.1142/s0218127400000840
https://doi.org/10.1142/s0218127400000840
https://doi.org/10.1007/978-1-4612-1828-9
https://doi.org/10.1007/978-1-4612-1828-9
https://doi.org/10.1113/jphysiol.1948.sp004260
https://doi.org/10.1113/jphysiol.1952.sp004764
https://doi.org/10.1007/978-0-387-87708-2
https://doi.org/10.1007/978-0-387-87708-2
https://doi.org/10.1162/neco.1996.8.5.979
https://doi.org/10.1016/j.camwa.2012.01.005
https://doi.org/10.1016/j.camwa.2012.01.005
https://doi.org/10.1038/nn.2212
https://doi.org/10.1016/j.apm.2019.12.008
https://doi.org/10.1016/j.apm.2019.12.008
https://doi.org/10.1016/j.cnsns.2013.10.032
https://doi.org/10.1007/s11571-013-9273-x
https://doi.org/10.1142/s0218127417502091
https://doi.org/10.1142/s0218127417502091
https://doi.org/10.1103/physrevlett.100.128103
https://doi.org/10.1016/j.mejo.2019.05.003
https://doi.org/10.1016/j.mejo.2019.05.003
https://doi.org/10.1038/s41598-019-52061-4
https://doi.org/10.1038/s41598-019-52061-4
https://doi.org/10.1002/(sici)1098-2426(199903)15:2<201::aid-num5>3.0.co;2-h


Eur. J. Math. Anal. 10.28924/ada/ma.3.2 27

[33] R.E. Mickens, A nonstandard finite difference scheme for a Fisher PDE having nonlinear diffusion, Computers Math.Appl. 45 (2003) 429–436. https://doi.org/10.1016/s0898-1221(03)80028-7.[34] R. Anguelov, J.M.S. Lubuma, Contributions to the mathematics of the nonstandard finite difference method andapplications, Numer. Meth. Part. Differ. Equ. 17 (2001) 518–543. https://doi.org/10.1002/num.1025.[35] C. Morris, H. Lecar, Voltage oscillations in the barnacle giant muscle fiber, Biophys. J. 35 (1981) 193–213. https:
//doi.org/10.1016/s0006-3495(81)84782-0.[36] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. 1 (1961)445–466. https://doi.org/10.1016/s0006-3495(61)86902-6.[37] T. Azizi, R. Mugabi, The phenomenon of neural bursting and spiking in neurons: Morris-Lecar model, Appl. Math.11 (2020) 203–226. https://doi.org/10.4236/am.2020.113017.[38] T. Azizi, B. Alali, Impact of Chloride Channel on Spiking Patterns of Morris-Lecar Model, Appl. Math. 11 (2020)650–669. https://doi.org/10.4236/am.2020.117044.[39] T. Azizi, Mathematical modeling with applications in biological systems, physiology, and neuroscience, KansasState University, 2021.[40] T. Azizi, B. Alali, G. Kerr, Mathematical modeling: with applications in physics, biology, chemistry, and engineering,Edition-2, (2021) 1–117. https://doi.org/10.9734/bpi/mono/978-93-91312-16-9.[41] R.K. Upadhyay, A. Mondal, Dynamics of fractional order modified Morris-Lecar neural model, Netw. Biol. 5 (2015)113-136.[42] G. Qi, Y. Wu, J. Hu, Abundant firing patterns in a memristive Morris–Lecar neuron model, Int. J. Bifurcation Chaos.31 (2021) 2150170. https://doi.org/10.1142/s0218127421501704.[43] M. Xing, X. Song, H. Wang, Z. Yang, Y. Chen, Frequency synchronization and excitabilities of two coupled hetero-geneous Morris-Lecar neurons, Chaos Solitons Fractals. 157 (2022) 111959. https://doi.org/10.1016/j.chaos.
2022.111959.[44] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations,J. Comput. Appl. Math. 172 (2004) 65–77. https://doi.org/10.1016/j.cam.2004.01.033.[45] D.M. Grobman, Homeomorphism of systems of differential equations, Dokl. Akad. Nauk SSSR 128 (1959) 880–881[46] P. Hartman, On local homeomorphisms of Euclidean spaces, Bol. Soc. Mat. Mexicana, 5 (1960), 220-241.[47] P. Hartman, A lemma in the theory of structural stability of differential equations, Proc. Amer. Math. Soc. 11 (1960)610–620. https://doi.org/10.1090/s0002-9939-1960-0121542-7.

https://doi.org/10.28924/ada/ma.3.2
https://doi.org/10.1016/s0898-1221(03)80028-7
https://doi.org/10.1002/num.1025
https://doi.org/10.1016/s0006-3495(81)84782-0
https://doi.org/10.1016/s0006-3495(81)84782-0
https://doi.org/10.1016/s0006-3495(61)86902-6
https://doi.org/10.4236/am.2020.113017
https://doi.org/10.4236/am.2020.117044
https://doi.org/10.9734/bpi/mono/978-93-91312-16-9
https://doi.org/10.1142/s0218127421501704
https://doi.org/10.1016/j.chaos.2022.111959
https://doi.org/10.1016/j.chaos.2022.111959
https://doi.org/10.1016/j.cam.2004.01.033
https://doi.org/10.1090/s0002-9939-1960-0121542-7

	1. Introduction
	2. Grünwald-Letinkov approximation
	3. Description of model equations
	3.1. Fractional Morris-Lecar model
	3.2. Local stability analysis of fractional order Morris-Lecar model

	4. Numerical results
	5. Discussion
	References

