The Phenomenon of Neural Bursting and Spiking in Neurons: Morris-Lecar Model

Tahmineh Azizi¹, Robert Mugabi²

¹Institute of Computational Comparative Medicine, Department of Anatomy and Physiology, Department of Mathematics, Kansas State University, Manhattan, KS, USA.
²Department of Diagnostic Medicine, and Pathobiology, College of Veterinary Medicine, Kansas State University, Manhattan, KS, USA.

DOI: 10.4236/am.2020.113017                

Abstract

In this study, we explore the interesting phenomenon of firing spikes and complex dynamics of Morris-Lecar model. We consider a set of parameters such that the model exhibits a wide range of phenomenons. We investigate the influences of injected current and temperature on the spiking dynamics of Morris-Lecar model. Moreover, we study bifurcations, and computational properties of this neuron model. Also, we define a bound (Max and Min voltage) for membrane potential and a certain voltage value or threshold for firing the spikes. Studying the two co-dimension bifurcations demonstrates much more complicated behaviors for this single neuron model. We also describe the phenomenon of neural bursting, and investigate the dynamics of Morris-Lecar model as a square-wave burster, elliptic burster and parabolic burster.

Keywords

Morris-Lecar Model, Neural Bursting, Firing Spike, Bifurcation

1. Introduction

During recent decades, understanding the brain function and exploring its molecular and cellular mechanisms were one of the greatest challenges in different fields of science. Historically, most of the researches in neuroscience focused on only neuronal circuits and synaptic organizations. Indeed, the neurons without considering their electrophysiological properties were divided into excitatory and inhibitory neurons, and sometimes they had been counted to be identical to those in Hodgkin-Huxley’s squid axon [1] [2] [3] [4]. In 1948 Hodgkin injected a dc-current of varying amplitude, and discovered that some preparations could show repetitive spiking activities with arbitrarily low frequencies, while the others discharged in a narrow frequency band [1] [5] [6] [7]. However, this finding was widely ignored by the scientists until 1989 when Rinzel and Ermentrout published a seminal paper and showed that the difference in behavior is because of different bifurcation mechanisms of excitability [1] [8] [9].

Non-linear dynamical system theory has a very important role in the computational neuroscience research [1] [2] [3] [4] [5] [10]. Izhikevich in [1] explains that how the transition in behavior of a neuron corresponds to a bifurcation from equilibrium to a limit cycle attractor. If we consider the injected current as a bifurcation parameter, when it is small, the cell remains quiescent. However, when the injected current increases, the cell switches to fire repetitive spikes [1] [2] [3] [4] [5] [10] [11]. In dynamical system theory, the qualitative change in the behavior of a system is called bifurcation. Indeed, when we change the amplitude of the bifurcation parameter (which in this case is the injected current), the cell undergoes a transition from quiescence to repetitive spiking. According to the type of bifurcation which happens for a neuron model, we can divide the neurons into different classes such as the class of excitability, and or we can discuss about the existence of threshold, all-or-none spikes, post-inhibitory rebound spikes, subthreshold oscillations, bistability of rest and spiking states [1]. For example, the neurons with supercritical and subcritical hopf bifurcations are called resonator and the neurons with saddle-node bifurcation or SNLC bifurcations are integrator [1].

In this paper, we study the interesting dynamics and fluctuations of spiking patterns of Morris-Lecar model which is a reduced version of Hodgkin-Huxley neuron model. For a certain range of parameters value, Morris-Lecar model exhibits different types of local bifurcations such as hopf bifurcation, saddle node on invariant limit cycles and homoclinic bifurcation. Moreover, we demonstrate a temperature bound and injected current range for spiking activity of the neuron. Also, we study co-dimension two bifurcations such as Bautin or generalized hopf and Bogdanov-Takens bifurcations and we present the normal form of these bifurcations as well. At the end, we look at the complicated dynamics of Morris-Lecar model as a burster.

2. Description of Model Equations

In 1981, Kathleen Morris and Harold Lecar introduced a simple model to generate the action potentials [12]. The Morris-Lecar model describes the electrical activities of neurons with a system of two non linear ordinary differential equations and includes different channels. This model is a reduction version of the four dimensional Hodgkin-Huxley model keeping the main properties of spike generations with much simpler mathematical and computational analysis [12] [13]. The Morris-Lecar model consists of three channels a potassium channel, a leak and a calcium channel and has the following form

$\{\begin{array}{l}{C}_M\frac{\text{d}V}{\text{d}t}=I_{app}-g_L\left(V-E_L\right)-g_{K}n\left(V-E_K\right)-g_{Ca}{m}_{\infty }\left(V\right)\left(V-E_{Ca}\right)\\ =I_{app}-I_{ion}\left(V,n\right),\\ \frac{\text{d}n}{\text{d}t}=\varphi \left(n_{\infty }\left(V\right)-n\right)/{\tau }_n\left(V\right),\end{array}$ (1)

where

$m_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_1\right)/V_2\right)\right],$ (2)

${\tau }_n\left(V\right)=1/\cosh \left(\left(V-V_3\right)/\left(2V_4\right)\right),$ (3)

$n_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_3\right)/V_4\right)\right].$ (4)

and

$I_{ion}\left(V\mathrm{,}n\right)=g_L\left(V-E_L\right)+g_{K}n\left(V-E_K\right)+g_{Ca}{m}_{\infty }\left(V\right)\left(V-E_{Ca}\right)$ (5)

where V demonstrates membrane potential, and n the activation variable of the persistent $K^{+}$ current, so it is a two-dimensional vector $\left(V\mathrm{,}n\right)$. $E_K$, $E_{Ca}$, and $E_L$ denote the Nernst equilibrium potentials. $I_{app}$ demonstrates the injected current and $I_{ion}$ the ionic current. Parameter $\varphi$ is a temperature factor. $g_L$ is leak membrane conductance, $g_K$ is potassium membrane conductance and $g_{Ca}$ is calcium membrane conductance. Moreover, $C_M$ is the total membrane capacitance. Also, the voltage-sensitive steady-state activation function $m_{\infty }\left(V\right)$ and $n_{\infty }\left(V\right)$, and the time constant ${\tau }_n\left(V\right)$ can be measured experimentally. A useful way to demonstrate the electrical properties of a neuron is using an equivalent circuit as we can see in Figure 1 [1]. In this case, the total current has the form

$I=C\stackrel{\dot{}}{V}+I_{Ca}+I_K$ Kirchhoff’s Law

where, $I_k=g_k\left(V-E_K\right)$, $I_{Ca}=g_{Ca}\left(V-E_{Ca}\right)$ are the major ionic currents. Therefore

$C\stackrel{\dot{}}{V}=I-I_{Ca}-I_K$

Also, $E_K<V<E_{Ca}$ where $I_{Ca}$ (inward current) is negative and also $I_K$ is positive. Basically, the inward currents depolarize the neuron and outward currents hyperpolarize it.

This simple model shows different types of dynamics such as hopf bifurcation, saddle node on invariant limit cycles and homoclinic bifurcation. In Table 1, we can see a list of parameters that cause three types of dynamics in Morris-Lecar model [10]. As we see, in Table 1, $I_{app}$ has not been included since we consider it as a bifurcation parameter for bifurcation diagrams.

3. The Hopf Case

For the case of hopf parameters, the model (1) has one equilibrium point which is intersection point of V nullcline and n nullcline. Indeed, it is not possible to find the explicit solution for any of the cases in Table 1. In Figure 2, we have

Figure 1. Equivalent Circuit for model (1). $E_K$, $E_{Ca}$, and $E_L$ the Nernst equilibrium potentials. $I_{app}$ the injected current, $g_L$ leak membrane conductance, $g_K$ potassium membrane conductance, $g_{Ca}$ calcium membrane conductance, $C_M$ the total membrane capacitance.

Figure 2. Nullclines of model (1) with $I_{app}=0,25,50,100$.

demonstrated the nullclines of system (1) for $I_{app}=0,25,50,100$ and it shows that there is only one solution for the hopf case.

In Figure 3, we can see different behaviors of the neuron from resting to

Figure 3. Time series of model (1) with $I_{app}=20$ left, and $I_{app}=40$ right.

Table 1. Morris-lecar parameters [10].

spiking (the stable constant solutions are corresponding to the resting state and spiking state shows the existence of periodic solutions).

Figure 4 displays the occurence of limit cycle corresponding to hopf bifurcation with increasing $I_{app}$.

As it has been exhibited in Figure 5, some trajectories come back to the stable fixed point or resting state after a big counter clockwise excursion but there are some other trajectories that return to the resting state without firing a spike. In neuroscience point of view, this kind of behavior is called threshold for firing spike. As a result, model (1) has an evident threshold for firing spikes and it has been obtained numerically equals to $V=-20\text{mv}$. Also, depends on the initial values of membrane potential, the size of the excursion is different.

Using continuation software Matcont and choosing $I_{app}$ as a bifurcation parameter, we could detect numerically hopf bifurcation points. For $I_{app}=93.857569$

Figure 4. Trajectories of model (1) with $I_{app}=87$ (up, left), $I_{app}=88.25$ (up, right), $I_{app}=88.3$ (down, left), $I_{app}=90$ (down, right).

Figure 5. Threshold for firing a spike in model (1).

and $\left(V\mathrm{,}n\right)=\left(-\mathrm{25.270122,0.139673}\right)$ we have the first hopf point with the first Lyapunov coefficient 5.220161, and for $I_{app}=212.018818$ and $\left(V\mathrm{,}n\right)=\left(\mathrm{7.800664,0.595491}\right)$ we can see the second hopf point with the first Lyapunov coefficient 5.451163. Here, the first Lyapunov coefficients for two hopf points are positive. Thus, there should exist an unstable limit cycle, bifurcating from the equilibrium and it indicates the appearance of subcritical hopf bifurcation. When the value for injected current $I_{app}$ is small we have a stable equilibrium point. When, we increase the value of injected current, the behavior of system changes and for $I_{app}=90$, we can see a limit cycle appears. This qualitative changing that causes producing a limit cycle attractor from a stable equilibrium point is called hopf bifurcation [14] [15]. In model (1) the equilibrium point is a stable focus that has a pair of complex conjugate eigenvalues with negative real part. With increasing the injected current, the real part of the eigenvalues changes from negative to zero and with further increasing, to positive. It means that the stable focus loses its stability and a limit cycle appears. With further increasing of the injected current, the amplitude of the limit cycle also increases.

For $\left(V\mathrm{,}n\right)=\left(\mathrm{7.800664,0.595491}\right)$ the eigenvalues are ${\lambda }_{\mathrm{1,2}}=4.65764\pm i\left(0.148602\right)$. Because the real part of eigenvalues are positive, it implies that equilibrium point is unstable and with further steps we obtain the second hopf point $\left(V\mathrm{,}n\right)=\left(-\mathrm{25.270122,0.139673}\right)$ with eigenvalues ${\lambda }_{\mathrm{1,2}}=-8.58989\pm i\left(0.0797799\right)$. Since, the real part is negative, it means that the equilibrium point is stable.

The topological normal form for hopf bifurcation has the form:

$\stackrel{\dot{}}{r}=\alpha r+a{r}^3$ (6)

$\stackrel{\dot{}}{\theta }={\omega }_0+\beta r^2$ (7)

Here, $\beta$ does not have any dynamical effect. The normal form for first hopf point $\left(V\mathrm{,}n\right)=\left(\mathrm{7.800664,0.595491}\right)$ has the form below:

$\stackrel{\dot{}}{r}=4.65764r+5.451163r^3$ (8)

$\stackrel{\dot{}}{\theta }=0.148602$ (9)

and also

$\stackrel{\dot{}}{r}=4.65764r+5.451163r^3$ (10)

$\stackrel{\dot{}}{\theta }=-0.148602$ (11)

In normal form (8), $\stackrel{\dot{}}{\theta }>0$. $\theta$ is the angle of oscillations which is positive and increasing because the frequency of damped or sustained oscillations around this point ${\omega }_0$, is positive. But for normal form (10), $\stackrel{\dot{}}{\theta }<0$ which means that the frequency of damped or sustained oscillations around this point ${\omega }_0$, is negative.

To analysis the normal form (8), we have

$r\left(4.65764+5.451163r^2\right)=0$

Therefore,

$r=0,4.65764+5.451163r^2=0$ (12)

Here, $r=0$ is an equilibrium and because for $r=0$,

$\frac{\text{d}}{\text{d}r}\left[4.65764r+5.451163r^3\right]=4.65764>0$, as a result, this equilibrium is unstable. The equation $4.65764+5.451163r^2=0$ does not give us any periodic solutions or oscillatory behaviors.

For other hopf point $\left(V\mathrm{,}n\right)=\left(-\mathrm{25.270122,0.139673}\right)$, the normal form can be written as

$\stackrel{\dot{}}{r}=-8.58989r+5.220161r^3$ (13)

$\stackrel{\dot{}}{\theta }=0.0797799$ (14)

Here, $\stackrel{\dot{}}{\theta }>0$ means that the frequency of damped or sustained oscillations around this point, ${\omega }_0$, is positive and increasing. But for the othe normal form:

$\stackrel{\dot{}}{r}=-8.58989r+5.220161r^3$ 15)

$\stackrel{\dot{}}{\theta }=0.0797799$ (16)

$\stackrel{\dot{}}{\theta }<0$ which implies that the frequency of damped or sustained oscillations around this point ${\omega }_0$, is negative and decreasing.

To analysis the normal form (13):

$r\left(-8.58989+5.220161r^2\right)=0$

Therefore

$r=0,-8.58989+5.220161r^2=0$ (17)

Here, $r=0$ is an equilibrium and because for $r=0$

$\frac{\text{d}}{\text{d}r}\left[-8.58989r+5.220161r^3\right]=-8.58989<0$

As a result, this equilibrium is stable. The equation $-8.58989+5.220161r^2=0$ gives us the unstable periodic solution with amplitude:

$r=\sqrt{\frac{8.58989}{5.220161}}$ (18)

4. The SNLC Case

Fold bifurcation of limit cycle or SNLC happens when with increasing the injected current two limit cycles, one stable that is associated to the stable node and another one unstable limit cycle that is associated to a saddle point close to each other, collide and at the bifurcation moment, we only have one limit cycle. With further increasing the injected current, this limit cycle also disappears. Figure 6 exhibits the nullclines of system (1) using SNLC parameters value in Table 1. As we see, with increasing the injected current, the numbers of equilibrium point change from 3 to only one fixed point.

Figure 7 demonstrates the trajectories of system (1) with the SNLC parameters in Table 1, and for different values for injected current $I_{app}$.

Figure 6. The Nullclines of model (1) in the case of SNLC bifurcation.

Figure 7. Trajectories of model (1) in the case of SNLC bifurcation for (up, left) $I_{app}=5$, (up.right) $I_{app}=30$, (down, left) $I_{app}=42$, (down, right) $I_{app}=100$.

When we do the continuation of equilibrium point of SNLC case, we detect the first hopf point for $I_{app}=97.646159$ and $\left(V\mathrm{,}n\right)=\left(\mathrm{8.334122,0.396190}\right)$ and first Lyapunov coefficient 5.317042 and two complex conjugate eigenvalues ${\lambda }_{\mathrm{1,2}}=\left(6.13675\right)\pm i\left(0.252728\right)$. Here, similar to the hopf case, hopf bifurcation is subcritical since the first Lyapunov coefficient is positive and its normal form has the following form

$\stackrel{\dot{}}{r}=6.13675r+5.317042r^3$ (19)

$\stackrel{\dot{}}{\theta }=-0.252748$ (20)

Since, $\stackrel{\dot{}}{\theta }<0$, the frequency of damped or sustained oscillations around this point ${\omega }_0$, is negative and decreasing.

However, the analysis of normal form is just limited to the first equation of normal form (19):

$r\left(6.13675+5.317042r^2\right)=0$

Therefore,

$r=0,6.13675+5.317042r^2=0$ (21)

Here, $r=0$ is an equilibrium and because for $r=0$, $\frac{\text{d}}{\text{d}r}\left[6.13675r+5.317042r^3\right]=6.13675>0$, this equilibrium point is unstable. The equation $6.13675+5.317042r^2=0$ does not give us any periodic solutions or oscillatory behaviors.

Also, the continuation of equilibrium point gives a limit point bifurcation for $I_{app}=-9.949039$ and at the point $\left(V\mathrm{,}n\right)=\left(-\mathrm{4.048524,0.136501}\right)$ with the normal form coefficient $a=4.772860$ and the eigenvalues $\left({\lambda }_1\mathrm{,}{\lambda }_2\right)=\left(-\mathrm{5.61265,0.359026}\right)$. The normal form for this bifurcation can be written as:

$\stackrel{\dot{}}{V}=a\pm V^2$ (22)

and for this case, we can write the following normal form:

$\stackrel{\dot{}}{V}=4.772860-V^2$ (23)

Thus, $V=\pm \sqrt{4.772860}$. Here, we find an equilibrium manifold which is the parabola $4.772860=V^2$ and gives the appearance of two equilibria. The same analysis can be done for other normal form and it gives the parabola $-4.772860=V^2$ but in this case, we have a singularity of the fold type.

The third point which has been detected by continuation is a Neutral saddle corresponding to ${\lambda }_1+{\lambda }_2=0$ for $I_{app}=36.639168$ at $\left(V\mathrm{,}n\right)=\left(-\mathrm{23.5341020.016555}\right)$ with eigenvalues $\left({\lambda }_1\mathrm{,}{\lambda }_2\right)=\left(-\mathrm{0.0792728,0.0792728}\right)$. Further continuation gives us another limit point bifurcation for $I_{app}=39.963153$ and at the point $\left(V\mathrm{,}n\right)=\left(-\mathrm{29.389788,0.008514}\right)$ for which a stable and an unstable limit cycle collide and create a non hyperbolic cycle. The real eigenvalues are $\left({\lambda }_1\mathrm{,}{\lambda }_2\right)=\left(-\mathrm{0.0990488,}-1.10511\right)$. For this fold bifurcation, the normal form would be

$\stackrel{\dot{}}{V}=-5.212474+V^2$ (24)

Thus, $V=\pm \sqrt{5.212474}$. Here, we have an equilibrium manifold which is the parabola $5.212474=V^2$ and it implies to the appearance of two equilibria. The same analysis can be done for the second normal form and we get the parabola $-5.212474=V^2$ which gives a singularity of the fold type.

5. The Homoclinic Case

Saddle-Homoclinic bifurcation happens when a saddle point and a limit cycle collide as we are increasing the control parameter. At the moment of bifurcation, we have a periodic orbit such that its period goes to infinity and finally, this periodic orbit disappears. The trajectories of system (1) for homoclinic case have been demonstrated in Figure 8.

Also, as we can see in Figure 9 the numbers of fixed points of model (1) with

Figure 8. Trajectories of model (1) in the case of Saddle-Homoclinic bifurcation for (up, left) $I_{app}=0$, (up.right) $I_{app}=40$, (down, left) $I_{app}=50$, (down, right) $I_{app}=70$.

Figure 9. Nullclines of model (1) with $I_{app}=0,25,50,100$.

increasing $I_{app}=0$ to $I_{app}=100$, reduce from 3 equilibrium point to one equilibrium point. Indeed, changing the numbers of fixed point means that a qualitative changes or bifurcation happens in the system.

When we do continuation of equilibrium points, we detect a hopf point for $I_{app}=36.316266$ and $\left(V\mathrm{,}n\right)=\left(\mathrm{4.410760,0.294770}\right)$ with the first Lyapunov coefficient 3.765575. For this hopf point the eigenvalues are complex conjugate: ${\lambda }_1=-1.3955\pm i\left(0.378861\right)$. Therefore, this is a subcritical hopf bifurcation. The normal form of this bifurcation would be

$\stackrel{\dot{}}{r}=-1.3955r+3.765575r^3$ (25)

$\stackrel{\dot{}}{\theta }=-0.378861$ (26)

Since, $\stackrel{\dot{}}{\theta }<0$, the frequency of damped or sustained oscillations around this point ${\omega }_0$, is negative and decreasing.

Analysis of normal form gives us

$r\left(-1.3955+3.765575r^2\right)=0$

Therefore

$r=0,-1.3955+3.765575r^2=0$ (27)

Here, $r=0$ is an equilibrium and because for $r=0$, $\frac{\text{d}}{\text{d}r}\left[-1.3955r+3.765575r^3\right]=-1.3955<0$, this equilibrium point is stable. The equation $-1.3955+3.765575r^2=0$ gives us periodic solution or oscillatory behaviors with amplitude:

$r=\sqrt{\frac{1.3955}{3.765575}}$ (28)

when we continue along the curve of equilibrium points, we detect a limit point bifurcation for $I_{app}=-9.949039$ and at the point $\left(V\mathrm{,}n\right)=\left(-\mathrm{4.048518,0.136501}\right)$, with the normal form coefficient $a=3.297636$. In this case, the eigenvalues are $\left({\lambda }_1,{\lambda }_2\right)=\left(-1.3742,0.178384\right)$. The normal form for this limit point bifurcation has the following form

$\stackrel{\dot{}}{V}=3.297636-V^2$ (29)

Consequently, $V=\pm \sqrt{3.297636}$. The equilibrium manifold would be the parabola $3.297636=V^2$. The same analysis can be done for other normal form and it gives us the parabola $-3.297636=V^2$, but in this case, we have a singularity of the fold type.

We can define the type of the saddle homoclinic bifurcation by looking at the sign of the sum of the eigenvalues which is called saddle quantity. If ${\lambda }_1+{\lambda }_2<0$, then the saddle homoclinic bifurcation is called supercritical which is corresponding to the appearance or disappearance of a stable limit cycle, and if ${\lambda }_1+{\lambda }_2>0$, we have the subcritical saddle homoclinic orbit bifurcation and it is corresponding to the appearance or disappearance of a unstable limit cycle. As a result, since here ${\lambda }_1+{\lambda }_2<0$, we have a supercritical saddle homoclinic bifurcation.

By further continuation the equilibrium curve, we obtain a limit point bifurcation for $I_{app}=39.963153$ and at the point $\left(V\mathrm{,}n\right)=\left(-\mathrm{29.389788,0.008514}\right)$, with normal form coefficient $a=-4.526064$. The eigenvalues are $\left({\lambda }_1\mathrm{,}{\lambda }_2\right)=\left(-\mathrm{0.391585,1.96656}\right)$. Also, the normal form would be

$\stackrel{\dot{}}{V}=-4.526064+V^2$ (30)

and we have $V=\pm \sqrt{4.526064}$. Therefore, we have an equilibrium manifold which is the parabola $4.526064=V^2$ and two equilibria appear. The same analysis gives the parabola $-4.526064=V^2$ but in this case, we have a singularity of the fold type.

Because here ${\lambda }_1+{\lambda }_2>0$, we have the subcritical saddle homoclinic orbit bifurcation. In neuroscience point of view, the saddle homoclinic bifurcation implies to the appearance or disappearance of spiking behavior.

6. Co-Dimension Two Bifurcations

In this section, we focus on co-dimension two bifurcations with $I_{app}$ and $\varphi$ as bifurcation parameters. The purpose is exploring the influences of temperature and injected current simultaneously on neuron model (1). At first, we discover Bautin or generalized hopf (GH) points for which the first Lyapunov coefficient vanishes. Then, we study another type of co-dimension two bifurcation which is called Bogdanov-Takens (BT) for which the system has an equilibrium with a double zero eigenvalue [14].

We start with the continuation of hopf curve that bifurcates from the BT point and continues to reach a Bautin point named GH. With further continuation, we can see another BT point after the second GH. Here, two GH points are non degenerate because the second Lyapunov coefficients $l_2$ are non zero, $l_2=-1.556360$.

Here, for $\left(V\mathrm{,}n\right)=\left(-\mathrm{11.785736,0.285152}\right)$, and parameters value $\left(\varphi \mathrm{,}{I}_{app}\right)=\left(\mathrm{0.306345,124.470639}\right)$, we have the first Bautin point. In generalized hopf (Bautin) bifurcation, the equilibrium has a pair of complex conjugate eigenvalues and also at generalized hopf point the first Lyapunov coefficient for the hopf bifurcation becomes zero. The bifurcation point separates branches of subcritical and supercritical hopf bifurcations. For the parameter values near bifurcation, the system demonstrates two limit cycles that collide and disappear via a saddle-node bifurcation. Basically, in Bautin bifurcation we have changing the type of bifurcation from subcritical to super critical hopf bifurcation. It means that the sign of the first coefficient Lyapunove changes from positive to negative. When the first Lyapunov coefficient becomes zero, the bifurcation becomes degenerate and the dynamics of system satisfies the following topological normal form [14]:

$\stackrel{\dot{}}{z}=\left(\lambda +i\omega \right)z+l_{1}z{\left|z\right|}^2+l_{2}z{\left|z\right|}^4$

Here, $z\in ℂ$ is a complex number, $l_1$ called the first Lyapunov coefficient and $l_2$ called the second Lyapunove coefficient, and $\lambda$ is the real part of eigenvalues and $\omega$ demonstrates the imaginary part of eigenvalues. At the moment of Bautin bifurcation $\lambda =l_1=0$ and $l_2\ne 0$. Likewise, when $l_2>0$, we have subcritical Bautin bifurcation and when $l_2<0$, we have supercritical Bautin bifurcation. To begin bifurcation analysis, at first when $\lambda =0$, we have a hopf bifurcation and depending on the sign of $l_1$ we have supercritical or subcritical hopf bifurcation. Also, when the first and the second Lyapunov coefficients have different sign, the solutions branch collide and disappear at the half parabola $l_1^2-4\lambda l_2=0$ and they undergo the fold limit cycle.

Here, for $\left(V\mathrm{,}n\right)=\left(-\mathrm{11.785736,0.285152}\right)$ and parameters value $\left(\varphi \mathrm{,}{I}_{app}\right)=\left(\mathrm{0.306345,124.470639}\right)$, we have a Bautin point with eigenvalues ${\lambda }_{\mathrm{1,2}}=-5.0307+i\left(0.156687\right)$ and the system satisfies the following normal form:

$\stackrel{\dot{}}{V}=\left(-5.0307+i\left(0.156687\right)\right)V+\left(-1.556360\right)V{\left|V\right|}^4$ (31)

Because, $l_2<0$, we have supercritical Bautin bifurcation. The other Bautin point happens for $\left(V\mathrm{,}n\right)=\left(\mathrm{2.472096,0.507868}\right)$, and $\left(\varphi \mathrm{,}{I}_{app}\right)=\left(\mathrm{0.253856,165.685695}\right)$ with eigenvalues ${\lambda }_{\mathrm{1,2}}=-2.67147\pm i\left(0.28612\right)$ and the system satisfies the following normal form

$\stackrel{\dot{}}{V}=\left(-2.67147+i\left(0.28612\right)\right)V+\left(-3.920527\right)V{\left|V\right|}^4$ (32)

Because, $l_2<0$, we have supercritical Bautin bifurcation.

Generalized hopf (Bautin) bifurcation in polar coordinates has the following normal form [14]:

$\stackrel{\dot{}}{r}=r\left(l_1+l_2{r}^2-r^4\right),$

In our simulation, the curve LPC corresponds to the saddle-node bifurcation of periodic orbits. As we can see in Figure 10, for $\varphi =0.30634507$, we have Limit point cycle with Normal form coefficient = 1.604795.

Moreover, from Figure 11, it can be easily observed that, Bogdanov-Takens bifurcation can be located along a hopf bifurcation curves, and as we approach to Bogdanov-Takens point, 2 purely imaginary eigenvalues collide and we have a double zero eigenvalue [14] [15].

Bogdanov-Takens bifurcation occurs when an equilibrium undergoes hopf bifurcation and saddle-node bifurcation simultaneously and also it occurs when we have at least a two-dimensional system. In this case, the Jacobian matrix of an equilibrium has these properties: $det\left(J\right)=0$ corresponding to saddle-node bifurcation, and $\text{tr}\left(J\right)=0$ corresponding to hopf bifurcation and it has the form:

${J|}_{\left(k\mathrm{,0}\right)}=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$

Because of these two conditions, Bogdanov-Takens is a codimension two bifurcation that has the following normal form

Figure 10. Continuation of equilibrium point in generalized hopf bifurcation with $\varphi =0.30634507$.

Figure 11. Continuation of equilibrium point in generalized hopf bifurcation with $\varphi =0.30634507$.

$\stackrel{\dot{}}{u}=v$

$\stackrel{\dot{}}{v}=a+bu+u^2+\sigma uv$

where, $a\mathrm{,}b$ are the normal form coefficients, and the parameter $\sigma$ takes the values 1 and −1, negative shows that it is supercritical and positive, when it is sub critical Bogdanov-Takens bifurcation. Two Bogdanov-Takens points in Figure 11 are: $\left(V\mathrm{,}n\right)=\left(-\mathrm{28.744348,0.114090}\right)$, $\left(\varphi \mathrm{,}{I}_{app}\right)=\left(\mathrm{0.000000,83.645532}\right)$ with $\left(a\mathrm{,}b\right)=\left(-\mathrm{5.341083,}-1.363867\right)$ and the second Bogdanov-Takens bifurcation happens at $\left(V\mathrm{,}n\right)=\left(\mathrm{8.717678,0.610127}\right)$, with the parameter values $\left(\varphi \mathrm{,}{I}_{app}\right)=\left(-\mathrm{0.000000,222.452534}\right)$ with $\left(a\mathrm{,}b\right)=\left(-\mathrm{1.185987,3.751062}\right)$. Finally, we compared the effect of injected current and temperature in Figure 12 and Figure 13. In Figure 12, we can easily find a lower bound and upper bound for injected current and a maximum and minimum voltage bound corresponding to spiking activity of this single neuron. Also, Figure 13 gives us a range for temperature and a maximum and minimum voltage bound corresponding to firing spike for Morris-Lecar model.

7. Bursting Behaviors of the Morris-Lecar Model

For some neurons that have spiking behavior, by applying some changes, they may also exhibit bursting behavior. For a neuron with ability to fire the spike, by adding a slow resonant current or gating variable we can change the neuron state to be a burster. The reason for this type of behavior is modulating the spiking and slow activity by the help of a slow negative feedback. Using the slow parameters, a burster can control the fast subsystem that has spiking state. Classification of bursters depends on the type of bifurcation of equilibrium points and limit cycles [10].

Figure 12. Continuation of Limit point cycles. Influence of injected current on neuron activities, maximum and minimum voltage bound.

Figure 13. Continuation of Limit point cycles. Influence of temperature on neuron activities, maximum and minimum voltage bound.

7.1. Morris-Lecar Model as a Square-Wave Burster

The first type of bursting is square wave bursting which has two important properties [10]:

1) The repetitive spikes at membrane potential is more depolarized than the silent state.

2) The frequency of spiking decreases during the spiking state.

Bursting occurs for systems with at least three dimension. For Morris-Lecar model, we consider $I_{app}$ decreases during the repetitive firing state process and increases during the silent state. Then, this burster demonstrates slow negative feedback together with hysteresis in the fast dynamics which specifically happens for square-wave bursting. For this case, we add a calcium dependent potassium current and the system obtains the form [10]:

$\{\begin{array}{l}{C}_M\frac{\text{d}V}{\text{d}t}=I_{app}-g_L\left(V-E_L\right)-g_{K}n\left(V-E_K\right)-g_{Ca}{m}_{\infty }\left(V\right)\left(V-E_{Ca}\right)-I_{KCa},\\ \frac{\text{d}n}{\text{d}t}=\varphi \left(n_{\infty }\left(V\right)-n\right)/{\tau }_n\left(V\right),\\ \frac{\text{d}\left[Ca\right]}{\text{d}t}=\epsilon \left(-\mu I_{Ca}-K_{Ca}\left[Ca\right]\right).\end{array}$ (33)

where

$m_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_1\right)/V_2\right)\right],$

${\tau }_n\left(V\right)=1/\cosh \left(\left(V-V_3\right)/\left(2V_4\right)\right),$

$n_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_3\right)/V_4\right)\right].$

where, $I_{KCa}$ demonstrates the calcium dependent potassium current and equals $I_{KCa}=g_{KCa}z\left(V-E_K\right)$. Here, $g_{KCa}$ is the maximal conductance for $I_{KCa}$ and z is a gating variable with a Hill-like dependence on the near membrane calcium concentration, $\left[Ca\right]$, and $z=\frac{{\left[Ca\right]}^p}{{\left[Ca\right]}^p+1}$. Without loss of generality, we assume $p=1$. The last equation of system (33) is a balance equation for $\left[Ca\right]$. The parameter $\mu$ has been used to convert current into a concentration flux and includes the ratio of the cell’s surface area to the calcium compartment’s volume. The parameter $K_{Ca}$ implies to the calcium removal rate and $\epsilon$ represents the ratio of free to total calcium in the cell. Because calcium is highly buffered, $\epsilon$ is small and the calcium dynamics is slow. The two first equations in of system (33) are called the fast subsystem and the third equation is called the slow equation [10].

Here, $I_{KCa}$ called outward current. If conductance $g_{KCa}z$ is large, the cell has hyper-polarization state which is corresponding to resting behavior. Conversely, if $g_{KCa}z$ is small, the cell fires spikes. We have demonstrated the model (33) as a circuit in Figure 14. Also, we have the required bursting parameters for different types of bursters in Table 2 [10]. Moreover, we have demonstrated the dynamics of this burster for different $I_{app}$ in Figure 15.

7.2. Morris-Lecar Model as an Elliptic burster

For Morris-Lecar model as an elliptic burster, we used the model (33) with parameters for Elliptic bursting as we have in Table 2. Also, we have demonstrated the dynamics of Morris-Lecar model as an elliptic burster for different $I_{app}$ in Figure 16.

7.3. Morris-Lecar Model as a Parabolic Burster

Unlike two previous bursters which we need only one slow variable for bursting

Figure 14. Equivalent Circuit for model (33). $E_K$, $E_{Ca}$, and $E_L$ the Nernst equilibrium potentials. $I_{app}$ the injected current, $g_L$ leak membrane conductance, $g_K$ potassium membrane conductance, $g_{Ca}$ calcium membrane conductance, $C_M$ the total membrane capacitance.

Figure 15. Continuation of equilibrium points for $I_{app}=0,50,100,150,200$ and occurrence of cusp bifurcation with considering 2 free parameters, the horizontal curve corresponding to co-dimension two bifurcation and the vertical curves are corresponding to co-dimension one bifurcation with increasing $I_{app}$ from left to right. Morris-Lecar model as a Square-Wave burster.

Figure 16. Continuation of equilibrium points for $I_{app}=-50,0,45,100,150,200,250,300$, the horizontal curve corresponding to co-dimension two bifurcation and the vertical curves are corresponding to co-dimension one bifurcation with increasing $I_{app}$ from left to right. Morris-Lecar model as an Elliptic burster.

Table 2. Bursting parameters [10].

behavior and the occurrence of the bistability in time series for fast subsystem, in parabolic burster, we need at least two slow variables and the bursting is not because of the bistability and hysteresis loop. In parabolic burster, the model has the form:

$\{\begin{array}{l}{C}_M\frac{\text{d}V}{\text{d}t}=I_{app}-g_L\left(V-E_L\right)-g_{K}n\left(V-E_K\right)-g_{Ca}{m}_{\infty }\left(V\right)\left(V-E_{Ca}\right)-I_{KCa}\\ \frac{\text{d}n}{\text{d}t}=\varphi \left(n_{\infty }\left(V\right)-n\right)/{\tau }_n\left(V\right)\\ \frac{\text{d}\left[Ca\right]}{\text{d}t}=\epsilon \left(-\mu I_{Ca}-K_{Ca}\left[Ca\right]\right)\\ \frac{\text{d}s}{\text{d}t}=\epsilon \left(s_{\infty }\left(V\right)-s\right)/{\tau }_s\end{array}$ (34)

where

$m_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_1\right)/V_2\right)\right]$

${\tau }_n\left(V\right)=1/\cosh \left(\left(V-V_3\right)/\left(2V_4\right)\right)$

$n_{\infty }\left(V\right)=\frac{1}{2}\left[1+\tanh \left(\left(V-V_3\right)/V_4\right)\right]$

$s_{\infty }\left(V\right)=0.5\left(1+\tanh \left(V-12\right)/24\right)$

and calcium dependent potassium current $I_{KCa}$ is $I_{KCa}=g_{KCa}z\left(V-E_K\right)$ and also a new calcium current $I_{Cas}s\left(V-E_{Ca}\right)$ that is depending on the gating variable s, and considering the parameters in Table 2 for parabolic bursting. Here, $V\mathrm{,}n$ are two fast variables and $\left[Ca\right]\mathrm{,}s$ are two slow variables. The circuit corresponding to this neuron is presented in Figure 17. Finally, we have demonstrated the dynamics of this burster for different $I_{app}$ in Figure 18.

Figure 17. Equivalent Circuit for model (34), Morris-Lecar model as a Parabolic burster. $E_K$, $E_{Ca}$, and $E_L$ the Nernst equilibrium potentials. $I_{app}$ the injected current, $g_L$ leak membrane conductance, $g_K$ potassium membrane conductance, $g_{Ca}$ calcium membrane conductance, $C_M$ the total membrane capacitance.

Figure 18. Continuation of equilibrium points for $I_{app}=\mathrm{0,50,100,150,200}$ and occurrence of the zero-hopf and Bogdanov-Takense and cusp bifurcation with considering 2 free parameters for model (34), the horizontal curve corresponding to co-dimension two bifurcation and the vertical curves are corresponding to co-dimension one bifurcation with increasing $I_{app}$ from left to right. Morris-Lecar model as a Parabolic burster.

8. Conclusion

Understanding of the structure of brain and its dynamics has been facilitated using computer simulations. During the recent decades, our understanding about brain dynamics and the mechanisms of different neuron cells has been greatly improved. Indeed, the field of computational neuroscience has been started with the work of Hodgkin and Huxley in 1952 using nonlinear partial differential equations. The Hodgkin-Huxley model and its reduction related models developed and improved many different areas in mathematics. Recently, dynamical systems theory and computational methods have been used frequently to study neuron activities in a many of neuronal models. The collaboration between experimentalists and theoreticians in analysis of neuronal models provides many progresses in the area of neuroscience [16] [17] [18]. In this paper we studied spiking dynamics of a single neuron model which is a reduction of well-known Hodgkin-Huxley model and consists of a system of ordinary differential equations. Depending on the different parameters value, the model reproduces quiescent, spiking and bursting activities. We numerically discovered the hopf bifurcation, SNLC bifurcation and homocinic bifurcation and we presented their normal form for each case separately. Through bifurcation analysis and continuation of equilibrium point, we explored the complicated dynamics which happened by changing the injected current or changing the temperature in this neuron model. We could find a range for spiking activities of the neuron for injected current $I_{app}$ and temperature $\varphi$. We also discovered co-dimension two bifurcations such as Bautin or generalized hopf, Bogdanov-Takens and limit point cycles and we demonstrated their normal forms. We also described the phenomenon of neural bursting, and we used the continuation method to discover different bifurcations for three types of bursting behaviors. We found two new co-dimension bifurcations compared to two-dimensional Morris-Lecar model, cusp bifurcation and zero-hopf bifurcation. We presented the Morris-Lecar bursting model for square-wave, elliptic and parabolic burster and we displayed the circuit model corresponding to each type of burster. Finding other types of bursting for Morris-Lecar model can be a new research project which needs to further study about the phenomenon of neuronal bursting.

Acknowledgements

This work was supported by the Institute of Computational Comparative Medicine (ICCM) and department of Mathematics of Kansas State University. With a special thanks to Dr. Majid Jaberi-Douraki for his full support.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Izhikevich, E.M. (2007) Dynamical Systems in Neuroscience. MIT Press, Cambridge. https://doi.org/10.7551/mitpress/2526.001.0001
[2]	Izhikevich, E.M. (2003) Simple Model of Spiking Neurons. IEEE Transactions on Neural Networks, 14, 1569-1572. https://doi.org/10.1109/TNN.2003.820440
[3]	Izhikevich, E.M. (2004) Which Model to Use for Cortical Spiking Neurons? IEEE Transactions on Neural Networks, 15, 1063-1070. https://doi.org/10.1109/TNN.2004.832719
[4]	Izhikevich, M. (2000) Neural Excitability, Spiking and Bursting. International Journal of Bifurcation and Chaos, 10, 1171-1266. https://doi.org/10.1142/S0218127400000840
[5]	Hoppensteadt, F.C. and Izhikevich, E.M. (2012) Weakly Connected Neural Networks. Nonlinear Analysis: Real World Applications, Vol. 126, Springer Science & Business Media, Berlin.
[6]	Hodgkin, A.L. (1948) The Local Electric Changes Associated with Repetitive Action in a non-Medullated Axon. The Journal of Physiology, 107, 165-181. https://doi.org/10.1113/jphysiol.1948.sp004260
[7]	Hodgkin, A.L. and Huxley, A.F. (1952) A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve. The Journal of Physiology, 117, 500-544. https://doi.org/10.1113/jphysiol.1952.sp004764
[8]	Rinzel, J. and Ermentrout, G.B. (1989) Analysis of Neural Excitability and Oscillations. MIT Press, Cambridge.
[9]	Rinzel, J. and Ermentrout, G.B. (1998) Analysis of Neural Excitability and Oscillations. Methods in Neuronal Modeling, 2, 251-292.
[10]	Ermentrout, G.B. and Terman, D.H. (2010) Mathematical Foundations of Neuroscience. Vol. 35, Springer Science & Business Media, Berlin. https://doi.org/10.1007/978-0-387-87708-2
[11]	Ermentrout, B. (1996) Type I Membranes, Phase Resetting Curves, and Synchrony. Neural Computation, 8, 979-1001. https://doi.org/10.1162/neco.1996.8.5.979
[12]	Morris, C. and Lecar, H. (1981) Voltage Oscillations in the Barnacle Giant Muscle Fiber. Biophysical Journal, 35, 193-213. https://doi.org/10.1016/S0006-3495(81)84782-0
[13]	FitzHugh, R. (1961) Impulses and Physiological States in Theoretical Models of Nerve Membrane. Biophysical Journal, 1, 445-466. https://doi.org/10.1016/S0006-3495(61)86902-6
[14]	Kuznetsov, Y.A. (2013) Elements of Applied Bifurcation Theory. Vol. 112, Springer Science & Business Media, Berlin.
[15]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. Vol. 2, Springer Science & Business Media, Berlin.
[16]	Gerstner, W. and Kistler, W.M. (2002) Spiking Neuron Models: Single Neurons, Populations, Plasticity. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511815706
[17]	Dayan, P., Abbott, L.F., et al. (2001) Theoretical Neuroscience. Vol. 806, MIT Press, Cambridge.
[18]	Nagumo, J., Arimoto, S. and Yoshizawa, S. (1962) An Active Pulse Transmission Line Simulating Nerve Axon. Proceedings of the IRE, 50, 2061-2070. https://doi.org/10.1109/JRPROC.1962.288235