∎

11institutetext: Leila Eftekhari 22institutetext: Department of Mathematics,
Tarbiat Modares University, Tehran, Iran
22email<EMAIL_ADDRESS>33institutetext: Soleiman Hosseinpour
44institutetext: Department of Applied Mathematics, Shahrood University of
Technology, Shahrood, Iran
44email<EMAIL_ADDRESS>55institutetext: Salvador Jiménez
66institutetext: Departamento de Matemática Aplicada a las TIC, Universidad
Politécnica de Madrid, Madrid, Spain
66email<EMAIL_ADDRESS>77institutetext: Moein Khalighi 88institutetext:
Department of Computing, University of Turku, Finland
88email<EMAIL_ADDRESS>
# Syntax Highlighting in LaTeX with the listings Package

writeLaTeX  Leila Eftekhari  Soleiman Hosseinpour  Moein Khalighi  Salvador
Jiménez

(Received: date / Accepted: date)

# Bifurcation analysis of a fractional-order Pinsky-Rinzel model

writeLaTeX  Leila Eftekhari  Soleiman Hosseinpour  Moein Khalighi  Salvador
Jiménez

(Received: date / Accepted: date)

###### Abstract

The present work describes a new fractional-order system of a two-compartment
CA3 hippocampal pyramidal cell, so-called the Pinsky-Rinzel model, with Caputo
fractional derivative. We investigate the transient of the solutions and the
general behavior of the system based on the bifurcation diagrams, in which
fractional derivative orders and currents injection are taken as bifurcation
parameters. Chaotic regions are obtained by tunning the values of the
fractional derivative orders and injection currents. The proposed numerical
approach helps to study the stability of the system under certain conditions.

###### Keywords:

CA3 hippocampal pyramidal neurons, Caputo fractional derivative, Bifurcation
analysis.

##  1 Introduction

The CA3 region in the brain system has received more attention in recent years
due to its particular role in memory encoding kandel2000principles . This area
is critical for initiating hippocampal interictal-like activity, which depends
on the presence of abundant recurrent excitatory connections between CA3
pyramidal cells miles1983single . CA3 pyramidal cells intrinsically display
various firing patterns at the cellular level, ranging from single action
potentials to complex bursts bib3 . Such bursting is important for place-cell
activity signal propagation and the induction of synaptic plasticity bib4 . On
the other hand, mathematical models aim to understand the behavior of CA3
pyramidal cells bib5 . In 1994, Pinsky and Rinzel developed a two-compartment
model for CA3 hippocampal pyramidal neurons in a guinea pig based on a complex
19 compartments Traub model bib1 . Their proposed model explains how
interactions between the somatic and dendritic compartments occur, thus
allowing a computational implementation and attracting many scholars’
attention bib4 .

The Pinsky-Rinzel model is a non-smooth, and classic methods’ mathematical
analysis of its dynamical behavior is inefficient. In 2001, Hahn and Drand
carried out a mathematical analysis of the dynamical properties of the model
based on Strogatz’s worksbib3 . They probed the changeover between resting,
bursting, and spiking states affected by increasing the amount of
extracellular potassium concentration. Due to the complexity and non-smooth
nature of the system, analyzing the dynamic scenery of the system was quite
difficult. For these reasons, L. A. Atherton and et al. modified the Pinsky-
Rinzel model bib3 . These models investigate both qualitative acontinuous
functionsnd quantitative analysis of the CA3 pyramidal cells, ranging from
detailed, multi-compartmental models down to single compartmental models. This
is important since it allows performing bifurcation analysis and stability
investigations and paves the way to implementing more numerical methods.

A bifurcation diagram is a good representation to identify the disappearance,
appearance, or stability changes like the limit sets (equilibrium points,
periodic or almost periodic solutions) as the bifurcation parameter varies.
This concept has been significantly studied in various fields during the last
decades due to its importance for understanding the behavior of dynamic
systems PhysRevResearch.2.023281 ; article1 ; article2 . In the field of
Neuroscience, authors in bib3  reported a bifurcation analysis of a two-
compartment hippocampal pyramidal cell model, which is the only bifurcation
analysis of the CA3 pyramidal cell model until now. They used a smooth system
to investigate the qualitative behavior of the original model. However, in the
current work, we generalized a CA3 pyramidal cell model, known as the Pinsky-
Rinzel model, with the concept of fractional calculus to study the behavior of
the action optional and its impact on the gating variables in the proposed
model. Fractional calculus is a generalization of ordinary derivatives and
integrals to non-integer order. It has been shown that fractional calculus
operators hold non-local properties and can represent memory effects. In
dynamical systems theory, these effects are described by explicit non-local
terms of past system states nagy2014efficient .

Many new fractional-order systems have been suggested in recent years that
they have rich and complex dynamics, such as periodic solutions, chaos
windows, and bifurcation. For instance, the Chua system with a derivative
order 2.7 makes chaos motion hartley1995chaos . Chaotic dynamics of a damped
van der Pol equation with fractional order are reported in chen2008chaotic .
Authors in gao2005chaos  investigated the nonlinear dynamics and chaos in a
fractional-order financial system. A periodically forced complex Duffing’s
oscillator was proposed, and chaos for the system was introduced in detail
chen2008nonlinear . Bifurcation, chaos control, and synchronization were
extended for a fractional-order Lorenz system with complex variables
gao2005chaos . In fractional-order systems, a bifurcation implies a
qualitative or topological change in dynamics with a variation of a system
parameter or derivative order, and bifurcation analysis becomes harder due to
the non-local property of the operator of fractional calculus. Concerning the
importance of bifurcation analysis in fractional-order systems, the number of
studies on fractional-order models is still limited in Neuroscience.

In this paper, motivated by the above discussion, we survey the solutions
(membrane potentials and currents) of the fractional-order Pinsky-Rinzel
model. We study the bifurcations of the model as the fractional derivative
order changes, and we identify regular and chaotic regimes. Next, we consider
the somatic and dendritic injections, ISapp, and IDapp , as bifurcation
parameters and explore the chaotic behavior of the system in each case.
Finally, we investigate the system’s stability condition and try to find some
conditions that make the equilibria asymptotically stable.

Figure 1: Schematic of the model and its dynamical behavior. (a) A two-
compartment model (LABEL:ca3) displaying currents from and to soma and
dendrite. The somatic compartment is combined with fast Sodium
$\mathrm{I}_{Na}$ and delayed rectifier Potassium $\mathrm{I}_{K_{DR}}$ and
leak current. The dendritic compartment has a persistent Calcium
$\mathrm{I}_{Ca}$, Calcium activated Potassium $\mathrm{I}_{K_{Ca}}$,
hyperpolarisation Potassium current $\mathrm{I}_{K_{AHP}}$, synaptic current
$I_{syn}$, and leak current. Leak currents have been omitted for clarity.
Electronic coupling between the two compartments is modeled using two
parameters; $\mathrm{g}_{c}$ is the strength of coupling, and $p$ is the
percentage of the total area in the somatic-like compartment.
$\mathrm{I}_{S_{app}}$ and $\mathrm{I}_{D_{app}}$ are currents injected
through an electrode to couple the currents between the two compartments. (b)
Non-smooth curve and smooth approximation of $\mathrm{V}_{\mathrm{s}}$ and
$\mathrm{V}_{\mathrm{d}}$ from the model (LABEL:ca3) for two particular values
of $$\mathrm{I}_{S_{app}}$=2.5$ and 0.75 when $$\mathrm{I}_{D_{app}}$=0$. When
the membrane potential is close to the resting potential of the cell, these
channels stay in a closed state, but when the membrane potential increases to
a certain threshold value, these channels start to open quickly. At this
moment, the action potential is fired. When the action potential dies, the
membrane enters a hyperpolarized state; during this stage, no action potential
can fire. This period is considered the refractory period.

##  2 Model

The following sections describe the main features of the CA3 hippocampal
pyramidal neurons modeled by an ordinary differential equation and show how to
incorporate non-integer order derivatives into this model.

###  2.1 Pinsky-Rinzel model

The Pinsky-Rinzel model is based on two compartments illustrated in Fig. 1a.
The CA3 mathematical model proposed in bib3 ; bib4 ; bib5  is described by an
integer-order system of differential equations

$\displaystyle$\mathrm{C}_{\mathrm{m}}$$\mathrm{V}_{\mathrm{s}}^{{}^{\prime}}$=-$\mathrm{I}_{Leak}$-$\mathrm{I}_{Na}$-$\mathrm{I}_{K_{DR}}$-\dfrac{$\mathrm{I}_{DS}$}{$\mathrm{p}$}+\dfrac{$\mathrm{I}_{S_{app}}$}{$\mathrm{p}$},$
(1)
$\displaystyle$\mathrm{C}_{\mathrm{m}}$$\mathrm{V}_{\mathrm{d}}^{{}^{\prime}}$=-$\mathrm{I}_{Leak}$-$\mathrm{I}_{Ca}$-$\mathrm{I}_{K_{Ca}}$-$\mathrm{I}_{K_{AHP}}$+\dfrac{$\mathrm{I}_{SD}$}{(1-$\mathrm{p}$)}+\dfrac{$\mathrm{I}_{D_{app}}$}{(1-$\mathrm{p}$)}.$

The currents of the model are functions of potentials as follows

$\mathrm{I}_{Na}$
$\displaystyle=$\mathrm{g}_{Na}${m_{\infty}^{2}}($\mathrm{V}_{\mathrm{s}}$)h($\mathrm{V}_{\mathrm{s}}$-$\mathrm{V}_{Na}$),$
$\mathrm{I}_{K_{DR}}$
$\displaystyle=$\mathrm{g}_{K_{DR}}$n($\mathrm{V}_{\mathrm{s}}$-$\mathrm{V}_{K}$),$
$\mathrm{I}_{Ca}$
$\displaystyle=$\mathrm{g}_{Ca}\mathrm{s}$^{2}($\mathrm{V}_{\mathrm{d}}$-$\mathrm{V}_{N}$),$
$\mathrm{I}_{K_{Ca}}$
$\displaystyle=$\mathrm{g}_{k_{Ca}}$C\leavevmode\nobreak\
\chi($\mathrm{C}$a)($\mathrm{V}_{\mathrm{d}}$-$\mathrm{V}_{Ca}$),$
$\mathrm{I}_{K_{AHP}}$
$\displaystyle=$\mathrm{g}_{K_{AHP}}$q\leavevmode\nobreak\
($\mathrm{V}_{\mathrm{d}}$-$\mathrm{V}_{K}$),$ $\mathrm{I}_{SD}$
$\displaystyle=-$\mathrm{I}_{DS}$=$\mathrm{g}_{c}$($\mathrm{V}_{\mathrm{d}}$-$\mathrm{V}_{s}$),$
$\mathrm{I}_{Leak}$
$\displaystyle=$\mathrm{g}_{L}$($\mathrm{V}$-$\mathrm{V}_{L}$).$

The activation and inactivation variables should satisfy these equations

$\displaystyle\omega^{{}^{\prime}}($\mathrm{V}$)$
$\displaystyle=\dfrac{\omega_{\infty}($\mathrm{V}$)-\omega}{\tau_{\omega}($\mathrm{V}$)},$
(2) $\displaystyle\omega_{\infty}($\mathrm{V}$)$
$\displaystyle=\dfrac{\alpha_{\omega}($\mathrm{V}$)}{\alpha_{\omega}($\mathrm{V}$)+\beta_{\omega}($\mathrm{V}$)},$
(3) $\displaystyle\tau_{\omega}($\mathrm{V}$)$
$\displaystyle=\dfrac{1}{\alpha_{\omega}($\mathrm{V}$)+\beta_{\omega}($\mathrm{V}$)},$
(4)

where, independly, we consider $\omega=h,n,s,m,C$ and $q$. The rate functions
are defined as follows

$\displaystyle\alpha_{m}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=\dfrac{0.32(-46.9-$\mathrm{V}_{\mathrm{s}}$)}{\exp(\frac{-46.9-$\mathrm{V}_{\mathrm{s}}$}{4})-1},$
$\displaystyle\beta_{m}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=\dfrac{0.28($\mathrm{V}_{\mathrm{s}}$+19.9)}{\exp(\frac{$\mathrm{V}_{\mathrm{s}}$+19.9}{5})-1},$
$\displaystyle\alpha_{n}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=\dfrac{0.016(-24.9-$\mathrm{V}_{\mathrm{s}}$)}{\exp(\frac{-24.9-$\mathrm{V}_{\mathrm{s}}$}{5})-1},$
$\displaystyle\beta_{n}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=0.25\exp(-1-0.025$\mathrm{V}_{\mathrm{s}}$),$
$\displaystyle\alpha_{h}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=0.128\exp\Big{(}\dfrac{-43-$\mathrm{V}_{\mathrm{s}}$}{18}\Big{)},$
$\displaystyle\beta_{h}($\mathrm{V}_{\mathrm{s}}$)$
$\displaystyle=\dfrac{4}{1+\exp(\frac{(-20-$\mathrm{V}_{\mathrm{s}}$)}{5})},$
$\displaystyle\alpha_{s}($\mathrm{V}_{\mathrm{d}}$)$
$\displaystyle=\dfrac{1.6}{1+\exp(-0.072($\mathrm{V}_{\mathrm{d}}$-5))},$
$\displaystyle\beta_{s}($\mathrm{V}_{\mathrm{d}}$)$
$\displaystyle=\dfrac{0.02($\mathrm{V}_{\mathrm{d}}$+8.9)}{\exp(\frac{($\mathrm{V}_{\mathrm{d}}$+8.9)}{5})-1},$
$\displaystyle\alpha_{C}($\mathrm{V}_{\mathrm{d}}$)$
$\displaystyle=\dfrac{(1-H($\mathrm{V}_{\mathrm{d}}$+10))\exp(\frac{($\mathrm{V}_{\mathrm{d}}$+50)}{11}-\frac{($\mathrm{V}_{\mathrm{d}}$+53.5)}{27})}{18.975}+H($\mathrm{V}_{\mathrm{d}}$+10)(2\exp(\frac{(-53.5-$\mathrm{V}_{\mathrm{d}}$)}{27})),$
$\displaystyle\beta_{C}($\mathrm{V}_{\mathrm{d}}$)$
$\displaystyle=(1-H($\mathrm{V}_{\mathrm{d}}$+10))(2\exp(\frac{(-53.5-$\mathrm{V}_{\mathrm{d}}$)}{27})-\alpha_{c}($\mathrm{V}_{\mathrm{d}}$)),$
$\displaystyle\alpha_{q}($\mathrm{C}\mathrm{a}$)$
$\displaystyle=\min(0.00002$\mathrm{C}\mathrm{a}$,0.01),$
$\displaystyle\beta_{q}($\mathrm{C}\mathrm{a}$)$ $\displaystyle=0.001,$
$\displaystyle{\chi}($\mathrm{C}\mathrm{a}$)$
$\displaystyle=\min\Big{(}\dfrac{$\mathrm{C}\mathrm{a}$}{250},1\Big{)},$

where $w=h,s,n,m$ are defined simply by continuous rate functions, while the
$C$, $q$, and $\chi$ are formulated as discontinuous rate functions where
$H(.)$ is the Heaviside step function. However, this is a non-smooth
approximation, and for bifurcation analysis, these discontinuous functions can
be converted to a smooth form bib3 :

$\displaystyle C_{\infty}($\mathrm{V}_{\mathrm{d}}$)=$
$\displaystyle\left(\dfrac{1}{1+\exp(\dfrac{-10.1-V_{d}}{0.1016})}\right)^{0.00925},$
$\displaystyle\tau_{C}($\mathrm{V}_{\mathrm{d}}$)=$ $\displaystyle
3.627\exp(0.03704V_{d}),$ $\displaystyle q_{\infty}($\mathrm{C}\mathrm{a}$)=$
$\displaystyle
0.7894\exp(0.0002726$\mathrm{C}\mathrm{a}$)-0.7292\leavevmode\nobreak\
\exp(-0.01672$\mathrm{C}\mathrm{a}$),$
$\displaystyle\tau_{q}($\mathrm{C}\mathrm{a}$)=$ $\displaystyle
657.9\exp(-0.02023$\mathrm{C}\mathrm{a}$)+301.8\exp(-0.002381$\mathrm{C}\mathrm{a}$),$
$\displaystyle\chi($\mathrm{C}\mathrm{a}$)=$ $\displaystyle
1.073\sin(0.003453$\mathrm{C}\mathrm{a}$+0.08095)+0.08408\sin(0.01634$\mathrm{C}\mathrm{a}$-2.34)$
$\displaystyle+$ $\displaystyle 0.01811\sin(0.0348\leavevmode\nobreak\
$\mathrm{C}\mathrm{a}$-0.9918),$

also we have

$\displaystyle\dfrac{$\mathrm{d}\mathrm{C}\mathrm{a}$}{$\mathrm{d}\mathrm{t}$}=-0.13$\mathrm{I}$_{$\mathrm{C}\mathrm{a}$}-0.075$\mathrm{C}\mathrm{a}$.$

The first minus sign is based on the convention that inward currents are
negative. The non-smooth and the smooth approximation demonstrate very similar
dynamic behaviors (Fig. 1b). In Appendix, the notations and maximal
conductance parameters are given in Table 1.

###  2.2 Fractional-order formulation

In this section, a time-dependent Pinsky-Rinzel model is introduced by using
fractional calculus. There are several definitions for fractional derivatives
and integrals for different purposes de2014review . The Caputo fractional
derivatives are more applicable due to providing initial conditions with clear
physical interpretation. It motivates us to propose the CA3 model in the sense
of the Caputo fractional time derivative, which holds non-local properties.

The Caputo fractional time derivative of order $\alpha$ is defined by bib8

$\displaystyle{}^{C}\\!D^{\alpha}f(t)=\begin{cases}I^{n-\alpha}D^{n}f(t),&\leavevmode\nobreak\
\leavevmode\nobreak\ n-1<\alpha<n,\\\ D^{n}f(t),&\leavevmode\nobreak\
\leavevmode\nobreak\ \alpha=n.\end{cases}$ (5)

where $n\in\mathbb{N}$, $D^{n}$ is the integer derivative of order $n$, and
$I^{n-\alpha}$ is the factional integral of order $(n-\alpha)$ over function
$f(t)$ that is defined by

$\displaystyle
I^{n-\alpha}f(t)=\dfrac{1}{\Gamma(n-\alpha)}\int_{0}^{t}f(\tau)(t-\tau)^{n-1-\alpha}d\tau.$
(6)

where $\Gamma(.)$ is Euler’s Gamma function and $f(t)\in C^{n}$. Replacing the
integer order derivatives of (1) by non-integer order derivatives causes
dimensional inconsistencies in the equation. Hence, we should insert some
parameters in the fractional order model to make the same unit of measure for
both sides of the equations. Westerlund and Ekstam bib26  have shown that
Jacques Curie’s empirical law for the current through capacitors and
dielectrics causes a capacitive current-voltage relationship for a non-ideal
capacitor

$I_{C}^{\alpha}=C_{m}(\alpha)\,{}^{C}\\!D^{\alpha}$

where $\alpha\in(0,1)$, $C_{m}(\alpha)$ indicates the fractional-order
capacitance with units
$\mathrm{(}\mathrm{A}\mathrm{/}\mathrm{V}\mathrm{)}\mathrm{s}^{\alpha}$ and
$\,{}^{C}\\!D^{\alpha}$ is the Caputo fractional-order differential operator.
Here, we generalize the Pinsky-Rinzel model (LABEL:ca3) by considering
fractional-order derivatives of the somatic ($\mathrm{V}_{\mathrm{s}}$) and
dendritic ($\mathrm{V}_{\mathrm{d}}$) membrance potentials

$\displaystyle
C_{m}(\alpha)\,{}^{C}\\!D^{\alpha}V_{s}=-$\mathrm{I}_{Leak}$-$\mathrm{I}_{Na}$-$\mathrm{I}_{K_{DR}}$-\dfrac{$\mathrm{I}_{DS}$}{$\mathrm{p}$}+\dfrac{$\mathrm{I}_{S_{app}}$}{$\mathrm{p}$},$
(7) $\displaystyle
C_{m}(\alpha)\,{}^{C}\\!D^{\alpha}V_{d}=-$\mathrm{I}_{Leak}$-$\mathrm{I}_{Ca}$-$\mathrm{I}_{K_{Ca}}$-$\mathrm{I}_{K_{AHP}}$+\dfrac{$\mathrm{I}_{SD}$}{(1-$\mathrm{p}$)}+\dfrac{$\mathrm{I}_{D_{app}}$}{(1-$\mathrm{p}$)}.$

where $C_{m}(\alpha)=\dfrac{\tau^{\alpha}}{R_{m}}$, $R_{m}$ is the membrane
resistance, and $\tau$ is the time constant bib27 . Also we rewrite

$\displaystyle\,{}^{C}\\!D^{\alpha}\omega($\mathrm{V}$)$
$\displaystyle=\dfrac{\omega_{\infty}($\mathrm{V}$)-\omega}{\tau_{\omega}($\mathrm{V}$)},$
(8)

where $\omega={h,s,n,m}$, $C$, $q$ and ${\chi}$ are defined simply by
continuous rate functions in the same way.

##  3 Results

Figure 2:  Time series plots of $\mathrm{V}_{\mathrm{s}}$ and
$\mathrm{V}_{\mathrm{d}}$ from model (LABEL:ca3) for two particular values of
$$\mathrm{I}_{S_{app}}$=2.5$ and 0.75 when $$\mathrm{I}_{D_{app}}$=0$ and time
= 500.

This section discusses some dynamical properties, transient, and possible
chaotic behavior, that the fractional-order Pinsky-Rinzel reduced model
modifies. Our numerical method is based on the predictor-corrector Adams-
Bashforth-Moulton method for fractional differential equations dithelm ; bib28
.

The numerical solutions of the model is provided for two types of interest
somatic current injection $$\mathrm{I}_{S_{app}}$=2.5$ and 0.75
$($\mu\mathrm{A}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$)$ (Fig.S3, S4) Results
display somatic membrane potential in two distinctive firing modes. The higher
currents produce regular spiking (Fig. 2a), whilst low somatic current
injections cause low-frequency bursts in the system (Fig. 2b). In the first
mode, the dynamic of variable $q$ governing activation of
$\mathrm{I}_{K_{AHP}}$ and calcium concentration stays at a high level to
prevent a dendritic spike (Fig. S1). In contrast, in the second mode, a burst
is operated when the slow, hyperpolarizing current activation variable
decreases below a threshold, causing a calcium spike in the dendrites (Fig.
S2).

The Action potential is produced in neurons , and the movement of ions does
propagation through a special kind of voltage-gated ion channels in the
membrane bib30 . when the membrane potential is close to the resting potential
of the cell these channels stay in a closed state but when the membrane
potential increases to a certain threshold value, these channels start to open
quickly. At this moment, the action potential is fired. When the action
potential dies, the membrane enters a hyperpolarized state; ; no action
potential can fire during this stage. This period is considered as the
refractory period.

###  3.1 Studying transient

Plotting the variables versus time can only show whether or not the solutions
have reached the asymptotic regime. To comprehend of the proposed model and
for future development in bifurcation analysis, we study transients of the
system.

Figure 3:  Time series plots of $\mathrm{V}_{\mathrm{s}}$,
$\mathrm{V}_{\mathrm{d}}$, Ca, h, q, n, s and c from model (7) when
$$\mathrm{I}_{S_{app}}$=2.5$, $\alpha=0.95$ and time=4000.

Here, the transient is a regime in which the solution has a qualitative
behavior quite different than what is observed afterward. Usually, this
corresponds to an initial “short” time window after which the qualitative
behavior of the system is more uniform and it is in opposition to a ”steady
state”.

As can shown in (Fig. 3), with $$\mathrm{I}_{S_{app}}$=2.5$, variables
$\mathrm{V}_{\mathrm{s}}$, $h$, $n$ have short transient. On the other hand,
$\mathrm{V}_{\mathrm{d}}$, $s$, and Ca do not have short transient and have
drastic amplitude changes for times up to $\sim 67$ (the distance between the
time of the big peaks, which is an estimation). Finally, after $t=400$, the
asymptotic regime seems to be obtained. the asymptotic regime, follows the
transient; the new qualitative regime of the solutions observed for longer
times and that one may represent the behavior as t goes to infinity.

Figure 4:  Time serious plots of $\mathrm{V}_{\mathrm{s}}$,
$\mathrm{V}_{\mathrm{d}}$, Ca, h, q, n, s and c from model (7) when
$$\mathrm{I}_{S_{app}}$=0.75$, $\alpha=0.95$ and time=4000.

Variables $c$ and $q$ do not present the initial burst, but they have a
transient such that $t=400$, which means for these variables the behavior
after 400 or even 500 is quite different from the one before, and this change
is even much more evident for q (the average value grows and seem to reach a
plateau after 400 or 500). Furthermore, the system seems to reach the
asymptotic behavior after the transient, as can be seen if we extend the plot
to higher times (an asymptotic behavior can be chaotic or regular). Variables
$V_{s}$ and h do not show a change, $V_{d}$ has an initial peak almost at
double amplitude than the rest, and the six following peaks appear after
shorter times than all the other peaks that appear later. Something similar
occurs for variable n; although the amplitude of the peaks is consistent in
the graphic the times at which they occur are much shorter below 400 than
afterward. For variable s, the initial peak is very high, the amplitude
decreases and there is a change of behavior at about $t\approx 300$ and
variable Ca also changes the behavior at about $t\approx 500$ ( Fig.3).

|
---|---
|
Figure 5:  $\mathrm{V}_{\mathrm{d}}$ versus $\mathrm{V}_{\mathrm{s}}$ model (7) when $$\mathrm{I}_{S_{app}}$=2.5$ and $\alpha=0.95$ |
---|---
Figure 6:  Attractor from model (7) via $\mathrm{V}_{\mathrm{s}}$,
$\mathrm{V}_{\mathrm{d}}$, $\alpha=0.95$, case $$\mathrm{I}_{S_{app}}$=2.5$
(left) and case $$\mathrm{I}_{S_{app}}$=0.75$ (right) .

Furthermore, something similar appears in the case of
$$\mathrm{I}_{S_{app}}$=0.75$; the transient does not show a difference in the
amplitude but the periodical behavior. The transient for this case may be
$[0,300]$. After that, the solutions seem to have a periodic behavior with
period $\sim 67$ (see Fig. 4). For more investigation, the final time is
extended up to $t=10000$. The final time has extended drastically because a
time-span $\Delta t=1000$ is not large enough to reveal the behaviour of this
system since the data show that there is some intrinsic time-scale of the
order of $t_{s}\approx 30$; $\Delta t\approx 33\leavevmode\nobreak\ t_{s}$ is
a little short. To make sure, the time initially is extended to $t=2000$, then
to $t=3000$, and finally to $t=10000$ (indeed, in the end, $t\approx
333\leavevmode\nobreak\ t_{s}$). In order to simplify the study, some graphics
of the first two variables are represented, $V_{s}$ and $V_{d}$ (see Fig. 5).
There is a transient (not very long; $t\approx 400$) after which the solution
is confined to values of $V_{d}$ in $[-5,50]$ (for $t\in[0,10000]$) instead of
$[-10,80]$ (for $t\in[1000,2000]$), otherwise, the solution gives similar
results. Since the number of lines is different for the case,
$t\in[2000,3000]$, to check that this does not go to a regular solution for
longer times, the final time is extended to $t=10000$. Then, the behavior is
the same; the system is chaotic (at least up to this time), and the final
graphic for times $t\in[500,10000]$ shows what appears to be a strange
attractor (Fig. 6). If we do the same for the case, ${I_{S_{app}}}=0.75$, we
see that the solution is regular. If we remove the initial transient, we
obtain a simple, closed curve (Fig. 6).

###  3.2 Bifurcation analysis

Bifurcations significantly aim to understand the dynamics of fractional-order
systems bib24 ; bib25 ; bib12 , for instance, the Chen model bib24a  and the
Hodgkin-Huxley model bib15 . The most remarkable point about bifurcation
diagrams is that each branch looks like the whole diagram. The diagram
contains many such copies of itself. This self-similarity can be seen to
repeat itself at ever finer resolutions. Such behavior is characteristic of
geometric entities called fractals devaney1990chaos .

In this section, bifurcation analysis of system (7) is carried out to
investigate the dynamical behavior of the fractional-order Pinsky-Rinzel model
(7) in terms of $$\mathrm{I}_{S_{app}}$,$\mathrm{I}_{D_{app}}$$, and order of
fractional derivative $\alpha$. Qualitative changes in the system dynamics are
called bifurcations, and the parameter values at which they occur are called
bifurcation points. We consider
$$\mathrm{I}_{S_{app}}$=2.5,$\mathrm{I}_{D_{app}}$=0,$\mathrm{g}_{\mathrm{c}}$=2.1$
and $\alpha\in I=(0,1)$. The fractional-order $\alpha$ is considered as the
bifurcation parameter. The bifurcation diagram demonstrates the chaotic
behavior of the system.

Figure 7:  Bifurcation diagram for $\alpha\in[0.7,1]$ and
$($\mathrm{I}_{S_{app}}$=2.5$, $$\mathrm{I}_{D_{app}}$=0,$ and
$$\mathrm{g}_{\mathrm{c}}\mathrm{=}\mathrm{2}\,\mathrm{1}$)$.

According to Fig. 7, the system is not chaotic for a range of $\alpha$, from
0.7 to $\sim$0.93. There is probably some transient behavior of the
trajectories. Subsequently, we turn $\mathrm{I}_{S_{app}}$ as the bifurcation
parameter when $$\mathrm{I}_{D_{app}}$=0,$\mathrm{g}_{\mathrm{c}}$=2.1$. Fig.
8 shows that the order of derivatives $\alpha$ alter the chaotic behaviors of
the system. In fact, by increasing the $\mathrm{I}_{S_{app}}$, chaos stops
sooner in the systems with fewer order values.

|
---|---
|
|
Figure 8: Bifurcation diagram when the parameter $\mathrm{I}_{S_{app}}$ varied
with the different values of the order $\alpha$.

Moreover, this fact is true when we consider the $\mathrm{I}_{D_{app}}$ as the
bifurcation parameter (Fig. 9).

|
---|---
|
|
Figure 9: Bifurcation diagram when the parameter $\mathrm{I}_{D_{app}}$ varied
with the different values of the order $\alpha$.

In Figs. 8 and 9, two distinct regions are still visible, and we can clearly
see chaotic as well as non-chaotic regions. The points where the regions
change indicate fractional Hopf bifurcations (a loss of the equilibrium’s
stability). As seen in Figs. 8 and 9, the chaotic region starts to fade by
reducing the value of $\alpha$ from 1 to 0.8. In addition, Figs. 8, 9, and 7
do not show any self-similarity.

###  3.3 Stability analysis

In the following, we focus on the stability of the suggested fractional-order
systems under specific conditions. Finding an analytical solution for the
proposed model seems impossible in general. In the literature, investigating
the stability for such large dynamics has always been challenging. However,
for our analysis, we take into account a numerical scheme. We study the
stability of the system based on a theorem for the Jacobian matrix of the
system. Because of the non-linear nature of the system, we then obtain the
equilibrium points numerically.

###### Definition 3.1

The autonomous nonlinear differential system in the sense of Caputo is given
as follows

${}^{C}D^{\alpha}x(t)=f(x(t)),$ (9)

where $x(t)\in\mathbb{R}^{n}$, $f(x(t))\in\mathbb{R}^{n\times 1}$ is a
nonlinear function vector, $0<\alpha<1$.

###### Definition 3.2

The point $\varepsilon^{*}$ is an equilibrium point of system (9), if and only
if $f(\varepsilon^{*})=0$.

###### Definition 3.3

The equilibrium point $\varepsilon^{*}$ of (9) is said to be

(1) stable if $\forall\epsilon>0$, there exists a $\delta>0$ such that
$\|x(t)-\varepsilon^{*}\|<\epsilon$ holds for $\forall
x_{0}\in\\{z:\|z-\varepsilon^{*}\|<\delta\\}$ and $\forall t>0$;

(2) asymptotically stable if the equilibrium point is locally stable and
$\lim_{t\rightarrow+\infty}x(t)=\varepsilon^{*}$.

###### Theorem 3.1

Consider the non-linear fractional-order autonomous system (9). Let
$\varepsilon^{*}$ is an equilibrium point of the nonlinear system (9), then
equilibrium point $\varepsilon^{*}$ is asymptotically stable if eigenvalues of
Jacobian matrix $\lambda(J(\varepsilon^{*}))$ of the system satisfy
$\left|\mathrm{arg}(\mathrm{spec}(J(\varepsilon^{*})))\right|>\dfrac{\alpha\pi}{2}.$
Where $\mathrm{spec}(J(\varepsilon^{*}))$ is the set of all eigenvalues of
$J(\varepsilon^{*})$kaslik2012nonlinear .

Since we cannot establish the stability conditions for this system in a
general form, we attempt to find a numerical technique for certain values. One
can consider any other values based on their needs. The related codes for the
implementation of the methods can be seen in Appendix.

We use Matlab software and the $fsolve$ command to obtain the equilibrium
points of the system. Based on Theorem (3.1) we also compute the Jacobian
matrix of the system to study the stability of the equilibrium points. Deu to
chaotic behaviour of the system in the interval [0.93,1] we choose
$\alpha=0.95$ and try to find some domain of ${I_{S_{app}}}$ and
${I_{D_{app}}}$ where the system (7) would be asymptotically stable. Here, we
evaluate two cases. First, we consider $\alpha=0.95\in[0.93,1],$
$$\mathrm{I}_{S_{app}}$=0.75,\leavevmode\nobreak\ $\mathrm{I}_{D_{app}}$=0,$
and $g_{c}=2.1$. The equilibrium point of the system is

$\varepsilon^{*}_{1}\approx[1.0225,0.9267,0.9950,0.0088,0.0149,0.0117,0.0456,0.5355],$

and the eigenvalues of the Jacobian matrix are as

$\mathrm{spec}(J(\varepsilon^{*}_{1}))\approx[-3.3066,-2.4601,-1.0242,-0.3830,-0.3123,-0.0749,0.0007+0.0005i,0.0007-0.0005i],$

For the second case, we consider
$$\mathrm{I}_{S_{app}}$=2.5,\leavevmode\nobreak\ $\mathrm{I}_{D_{app}}$=0,$
and $g_{c}=2.1$. in the case2, the equilibrium point is given by

$\varepsilon^{*}_{2}\approx[2.6161,2.4916,0.9929,0.0239,0.0170,0.0132,0.0448,0.6876],$

and the eigenvalues are as

$\mathrm{spec}(J(\varepsilon^{*}_{2}))\approx[-3.2804,-2.3217,-0.9954,-0.3453,-0.2856,-0.0748,0.0039,-0.0000]$

In both cases, $\lvert arg(\lambda_{i})\rvert$ for $i=1,\cdots,8$, do not
satisfying the conditions of the Theorem (3.1) and then the equilibria
$\varepsilon^{*}_{1}$ and $\varepsilon^{*}_{2}$ are not asymptotically stable.

Now we try to find some domains of $\mathrm{I}_{S_{app}}$ and
$\mathrm{I}_{D_{app}}$ near the standard value of $$\mathrm{I}_{S_{app}}$=2.5$
and $$\mathrm{I}_{D_{app}}$=0$ to make the system (7) asymptotically stable.
At first, we consider $\mathrm{I}_{S_{app}}$ as a parameter that changes in
the interval $[-4,4]$. Using Remark(3.1) and an increment equal to $0.0001$ in
$\mathrm{I}_{S_{app}}$ we found the interval $[-1.2579,0.0268]$ where the
system does not show chaotic behaviour and the equilibria $\varepsilon^{*}$ is
asymptoticalliy stable (see Figs.10 ).

---
Figure 10:  Memberance potential for $\alpha=0.95$,
$$\mathrm{I}_{D_{app}}$=0$, $$\mathrm{I}_{S_{app}}$=-1.2579$ (left),
$$\mathrm{I}_{S_{app}}$=0.0268$ (right).

Secondely, taking $\mathrm{I}_{D_{app}}$ as a parameter changing in the
interval $[-4,4]$ with $$\mathrm{I}_{S_{app}}$=2.5$ and
$$\mathrm{g}_{\mathrm{c}}$=2.1$, using the Remark(3.1) and an increment equal
to $0.0001$ in $\mathrm{I}_{D_{app}}$ we found the interval $[-4,-2.5471]$
where the system does not show chaotic behaviour and also the equilibria
$\varepsilon^{*}$ of the system is asmptoticalliy stable (see Fig.11).
Finally, we let $\mathrm{I}_{D_{app}}$ as parameter in the interval [-4, 4]
with $$\mathrm{I}_{S_{app}}$=0.75$ and $$\mathrm{g}_{\mathrm{c}}$=2.1$. Using
the Remark(3.1) and an increment equal to $0.0001$ in $\mathrm{I}_{D_{app}}$
we found the interval $[-4,-0.7449]$ where the system doesn’t show chaotic
behaviour and also the equilibria $\varepsilon^{*}$ is asmptoticalliy stable
(see Fig.11).

---
Figure 11:  Memberance potential for $\alpha=0.95$,
$$\mathrm{I}_{S_{app}}$=2.5$ and $$\mathrm{I}_{D_{app}}$=-2.5471$ (left),
$$\mathrm{I}_{S_{app}}$=0.75$ and $$\mathrm{I}_{D_{app}}$=-0.7449$ (right).

##  4 Discussion

This work investigated the bifurcation analysis of a fractional-order model of
a two-compartment CA3 hippocampal pyramidal cell. Chaotic regions were
achieved for different values of the fractional derivative order and different
injection currents. Since the membrane capacitance has been considered ideal,
the fractional-order model has deemed as a more accurate description of
physical processes underlying a long-range memory behavior. We can model the
dielectric losses in the membrane with fractional capacitance and observe a
change in the relative refractory period, as shown in Fig. S5. In fact, Fig.
S5 shows a comparison between fractional-order and integer-order cases. In
fractional case $\alpha=0.85$ the refractory period is smaller than the
integer case $\alpha=1$. Moreover, based on Figs. 8 and 9, it can be concluded
that the new fractional-order system’s chaotic region collapses when the
derivative order approaches from 1 to 0. The obtained results can be
considered as help to control relevant diseases caused by maximal injection
currents abnormality.

## Appendix

In this section, we list the notations and the value of the parameters used
for the numerical results (Table 1), and illustrate the supplementary figures
(Figs S1-S4).

##  5 Matlab code

⬇

clc

close all

% y=[Vs, Vd, h, n, s, c, q, Ca]

% Id=I_{D_{app}} and Is= I_{S_{app}}, Ild,Ils=I_{Leak}

y0= [-4.6 ;-4.5 ;0.999 ;0.001 ;0.009 ;0.007; 0.01; 0.2]

Id = 0:.01:4;

i=0

alpha=0.95

for ct = 1:numel(Id)

y(:,ct)=fsolve(@(y)myfun(y,Id(ct)),y0);

df=NumJacob(@(y)myfun(y,Id(ct)),y(:,ct));

EIG= eig(df);

angle(EIG);

idx = find(abs(angle(EIG))<=(alpha*pi/2));

if (numel(idx~=0))

disp([’is␣not␣asymptotically␣stable’]);

␣␣␣␣␣␣␣else

␣␣␣␣␣i=i+1

␣␣␣␣␣disp([’is␣asymptotically␣stable’])

␣␣end

end

⬇

function f=test_func(y,Id)

gL = .1; gNa = 30; gKDR =15; gCa = 10; gKAHP = 0.8; gKC =15;

VNa = 120; VCa = 140; VK = -15; VL = 0;

Isyn = 0; p = 0.5; alpha=.95; Rm=10; tm=30; Cm=tm^alpha/Rm;

Is=2.5; gc=2.1;

Ils = gL*(y(1)-VL);

minf = am(y(1))./(am(y(1))+bm(y(1)));

INa = gNa*minf.^2.*y(3).*(y(1)-VNa);

IKDR = gKDR*y(4)*(y(1)-VK);

Ild = gL*(y(2)-VL);

ICa = gCa*y(5).^2.*(y(2)-VCa);

K=(1.073*sin(0.003453*y(8)+0.08095)+ 0.08408*sin(0.01634*y(8)-2.34)+
0.01811*sin(0.0348*y(8)-0.9918));

IKC = gKC*y(6)*K*(y(2)-VK);

IKAHP = gKAHP*y(7).*(y(2)-VK);

f = [(-Ils - INa - IKDR + (gc/p)*(y(2)-y(1)) + Is/p)/Cm;…

(-Ild-ICa-IKAHP-IKC-Isyn/(1-p)+(gc/(1-p))*(y(1)-y(2))+Id/(1-p))/Cm;…

ah(y(1))*(1-y(3)) - bh(y(1))*y(3);…

an(y(1))*(1-y(4)) - bn(y(1))*y(4);…

as(y(2))*(1-y(5)) - bs(y(2))*y(5);…

(aq(y(8))-y(7))/bq(y(8));

-0.13*ICa - 0.075*y(8)];

return

function val = am(v)

v=v-60;

val = 0.32*(-46.9-v)/(exp((-46.9-v)/4)-1);

function val = bm(v)

v=v-60;

val = 0.28*(v+19.9)/(exp((v+19.9)/5)-1);

function val = an(v)

v=v-60;

val = 0.016*(-24.9-v)/(exp((-24.9-v)/5)-1);

function val = bn(v)

v=v-60;

val =0.25*exp(-1-0.025*v);

function val = ah(v)

v=v-60;

val = 0.128*exp((-43-v)/18);

function val = bh(v)

v=v-60;

val =4/(1+exp((-20-v)/5));

function val = as(v)

v=v-60;

val =1.6/(exp(-0.072*(v-5))+1);

function val = bs(v)

v=v-60;

val =0.02*(v+8.9)/(exp((v+8.9)/5)-1);

function val = aq(v)

val =(0.7894*exp(0.0002726*v))-(0.7292*exp(-0.01672*v));

function val = bq(v)

val =(657.9*exp(-0.02023*v))+(301.8*exp(-0.002381*v));

function val = ac(v)

v=v-60;

val =(1.0/(1+exp((-10.1-v)/0.1016)))^0.00925;

function val = bc(v)

v=v-60;

val =3.627*exp(0.03704*v);

Table 1: Notations of the model and maximal conductance parameters. Notations | Description | Values
---|---|---
$\mathrm{I}_{Leak}$ | Leak currents | —
$\mathrm{I}_{Na}$ | Sodium current | —
$\mathrm{I}_{K_{DR}}$ | Delayed rectifier K current | —
$\mathrm{I}_{DS}$ | Coupling current | —
$\mathrm{I}_{SD}$ | Coupling current | —
$\mathrm{I}_{S_{app}}$ | Somatic current injection | —
$\mathrm{I}_{Ca}$ | Dendritic calcium current | —
$\mathrm{I}_{K_{Ca}}$ | Voltage and calcium dependent potassium current | —
$\mathrm{I}_{K_{AHP}}$ | Ca dependent AHPK current | —
$\mathrm{I}_{D_{app}}$ | Dendritic current injection | —
$\mathrm{V}_{\mathrm{s}}$ | Somatic membrance potentials | —
$\mathrm{V}_{\mathrm{d}}$ | Dendritic membrance potentials | —
$\mathrm{V}_{Na}$ | Reversal potentials, relative to resting potential, for$Na^{+}$ | 60
$\mathrm{V}_{K}$ | Reversal potentials, relative to resting potential, for $K^{+}$ | -75
$\mathrm{V}_{Ca}$ | Reversal potentials, relative to resting potential, for $Ca^{+}$ | 80
$\mathrm{V}_{L}$ | Reversal potentials, relative to resting potential, related to $I_{L}$ | -60
$\mathrm{p}$ | The size of the axosomatic compartment as a proportion of the entire cell | 0.5
$\mathrm{1}\mathrm{-}\mathrm{p}$ | The size of the dendritic compartment as a proportion of the entire cell | 0.5
$\mathrm{C}_{\mathrm{m}}$ | The capacitance | 3
$\mathrm{g}_{Na}$ | Maximum conductance for Na | 30
$\mathrm{g}_{K_{DR}}$ | Maximum conductance for delayed rectifier K current | 15
$\mathrm{g}_{K_{Ca}}$ | Maximum conductance for voltage and Ca dependent K current (fast) | 15
$\mathrm{g}_{K_{AHP}}$ | Maximum conductance for Ca dependent AHPK current (slow) | 0.8
$\mathrm{g}_{Ca}$ | Maximum conductance for Ca | 10
$\mathrm{g}_{L}$ | Leakage conductance | 0.1
$\mathrm{g}_{c}$ | Strength of coupling | 2.1
Figure S1: Activation and inactivation variables when
$$\mathrm{I}_{S_{app}}$=2.5$ and $$\mathrm{I}_{D_{app}}$=0$. Figure S2:
Activation and inactivation variables when $$\mathrm{I}_{S_{app}}$=0.75$ and
$$\mathrm{I}_{D_{app}}$=0$. Figure S3: Time course of the currents in the
somatic compartment, from somatic to dendritic compartment, and in the
dendritic compartment when $$\mathrm{I}_{S_{app}}$=2.5$ and
$$\mathrm{I}_{D_{app}}$=0$.. Figure S4: Time course of the currents in the
somatic compartment, from somatic to dendritic compartment, and in the
dendritic compartment when $$\mathrm{I}_{S_{app}}$=0.75$ and
$$\mathrm{I}_{D_{app}}$=0$. Figure S5: refractory period ($t_{rp2}$) for
$\alpha=0.85$ (yellow) and refractory period ($t_{rp1}$) for $\alpha=1$ (red).

## References

  * [1] Laura A Atherton, Luke Y Prince, and Krasimira Tsaneva-Atanasova.  Bifurcation analysis of a two-compartment hippocampal pyramidal cell model.  Journal of computational neuroscience, 41(1):91–106, 2016.
  * [2] Donato Cafagna and Giuseppe Grassi.  Bifurcation and chaos in the fractional-order chen system via a time-domain approach.  International Journal of Bifurcation and Chaos, 18(07):1845–1863, 2008.
  * [3] Juhn-Horng Chen and Wei-Ching Chen.  Chaotic dynamics of the fractionally damped van der pol equation.  Chaos, Solitons & Fractals, 35(1):188–198, 2008.
  * [4] Wei-Ching Chen.  Nonlinear dynamics and chaos in a fractional-order financial system.  Chaos, Solitons & Fractals, 36(5):1305–1314, 2008.
  * [5] Marius-F. Danca and Nikolay Kuznetsov.  Matlab code for lyapunov exponents of fractional-order systems.  International Journal of Bifurcation and Chaos, 28(05):1850067, 2018\.
  * [6] Edmundo Capelas de Oliveira and José António Tenreiro Machado.  A review of definitions for fractional derivatives and integral.  Math. Probl. Eng., 2014, 2014.
  * [7] Robert L Devaney.  Chaos, fractals, and dynamics: computer experiments in mathematics.  Addison-Wesley Longman Publishing Co., Inc., 1990.
  * [8] Kai Diethelm and Neville J Ford.  Analysis of fractional differential equations.  Journal of Mathematical Analysis and Applications, 265(2):229–248, 2002.
  * [9] G Bard Ermentrout and David H Terman.  Mathematical foundations of neuroscience, volume 35.  Springer Science & Business Media, 2010.
  * [10] Gabbiani F and Cox S J.  Mathematics for Neuroscientists.  Elsevier, 2017.
  * [11] Xin Gao and Juebang Yu.  Chaos in the fractional order periodically forced complex duffing’s oscillators.  Chaos, Solitons & Fractals, 24(4):1097–1104, 2005.
  * [12] Roberto Garrappa.  Numerical solution of fractional differential equations: A survey and a software tutorial.  Mathematics, 6(2):16, 2018.
  * [13] Tom T Hartley, Carl F Lorenzo, and H Killory Qammer.  Chaos in a fractional order chua’s system.  IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 42(8):485–490, 1995.
  * [14] Podlobny I.  Fractional Differential Equations.  Academic press, 1999.
  * [15] Eric R Kandel, James H Schwartz, Thomas M Jessell, Steven Siegelbaum, A James Hudspeth, and Sarah Mack.  Principles of neural science, volume 4.  McGraw-hill New York, 2000.
  * [16] Eva Kaslik and Seenith Sivasundaram.  Nonlinear dynamics and chaos in fractional-order neural networks.  Neural Networks, 32:245–256, 2012.
  * [17] A. Khan.  Bifurcation analysis of a discrete-time two-species model.  Discrete Dynamics in Nature and Society, 2020:1–12, 01 2020.
  * [18] Xiaojun Liu, Ling Hong, Lixin Yang, and Dafeng Tang.  Bifurcations of a new fractional-order system with a one-scroll chaotic attractor.  Discrete Dynamics in Nature and Society, 2019:8341514, Jan 2019\.
  * [19] Richard Miles and Robert KS Wong.  Single neurones can initiate synchronized population discharge in the hippocampus.  Nature, 306(5941):371–373, 1983.
  * [20] AM Nagy and NH Sweilam.  An efficient method for solving fractional hodgkin–huxley model.  Physics Letters A, 378(30-31):1980–1984, 2014.
  * [21] Paul Redfern and Stephen Thesleff.  Action potential generation in denervated rat skeletal muscle: Ii. the action of tetrodotoxin.  Acta physiologica scandinavica, 82(1):70–78, 1971.
  * [22] M. Sayed, Abd allah Mousa, and Ibrahim Mustafa.  Stability and bifurcation analysis of a buckled beam via active control.  Applied Mathematical Modelling, 82, 02 2020.
  * [23] R. D. Traub, R. K. Wong, R. Miles, and H. Michelson.  A model of a ca3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances.  Journal of Neurophysiology, 66(2):635–650, 1991.  PMID: 1663538.
  * [24] Seth H Weinberg.  Membrane capacitive memory alters spiking in neurons described by the fractional-order hodgkin-huxley model.  PloS one, 10(5):e0126629, 2015.
  * [25] Svante Westerlund and Lars Ekstam.  Capacitor theory.  IEEE Transactions on Dielectrics and Electrical Insulation, 1(5):826–839, 1994.
  * [26] Can Xu, Xuebin Wang, and Per Sebastian Skardal.  Bifurcation analysis and structural stability of simplicial oscillator populations.  Phys. Rev. Research, 2:023281, Jun 2020.
  * [27] Yue Zhang, Kuanquan Wang, Yongfeng Yuan, Dong Sui, and Henggui Zhang.  Stability and bifurcation analysis of hodgkin-huxley model.  In 2013 IEEE International Conference on Bioinformatics and Biomedicine, pages 49–54. IEEE, 2013.
  * [28] Jan Čermák and Luděk Nechvátal.  Local bifurcations and chaos in the fractional rössler system.  International Journal of Bifurcation and Chaos, 28(08):1850098, 2018\.