ath/MathML"> θ ˙ < 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.220161{r}^{2}\right)=0$

Therefore

$r=0,\text{ }-8.58989+5.220161{r}^{2}=0$ (17)

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

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

As a result, this equilibrium is stable. The equation $-8.58989+5.220161{r}^{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,n\right)=\left(8.334122,0.396190\right)$ and first Lyapunov coefficient 5.317042 and two complex conjugate eigenvalues ${\lambda }_{1,2}=\left(6.13675\right)±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{˙}{r}=6.13675r+5.317042{r}^{3}$ (19)

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

Since, $\stackrel{˙}{\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.317042{r}^{2}\right)=0$

Therefore,

$r=0,\text{ }6.13675+5.317042{r}^{2}=0$ (21)

Here, $r=0$ is an equilibrium and because for $r=0$, $\frac{\text{d}}{\text{d}r}\left[6.13675r+5.317042{r}^{3}\right]=6.13675>0$, this equilibrium point is unstable. The equation $6.13675+5.317042{r}^{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,n\right)=\left(-4.048524,0.136501\right)$ with the normal form coefficient $a=4.772860$ and the eigenvalues $\left({\lambda }_{1},{\lambda }_{2}\right)=\left(-5.61265,0.359026\right)$. The normal form for this bifurcation can be written as:

$\stackrel{˙}{V}=a±{V}^{2}$ (22)

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

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

Thus, $V=±\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,n\right)=\left(-23.5341020.016555\right)$ with eigenvalues $\left({\lambda }_{1},{\lambda }_{2}\right)=\left(-0.0792728,0.0792728\right)$. Further continuation gives us another limit point bifurcation for ${I}_{app}=39.963153$ and at the point $\left(V,n\right)=\left(-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},{\lambda }_{2}\right)=\left(-0.0990488,-1.10511\right)$. For this fold bifurcation, the normal form would be

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

Thus, $V=±\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,n\right)=\left(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±i\left(0.378861\right)$. Therefore, this is a subcritical hopf bifurcation. The normal form of this bifurcation would be

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

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

Since, $\stackrel{˙}{\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.765575{r}^{2}\right)=0$

Therefore

$r=0,\text{ }-1.3955+3.765575{r}^{2}=0$ (27)

Here, $r=0$ is an equilibrium and because for $r=0$, $\frac{\text{d}}{\text{d}r}\left[-1.3955r+3.765575{r}^{3}\right]=-1.3955<0$, this equilibrium point is stable. The equation $-1.3955+3.765575{r}^{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,n\right)=\left(-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{˙}{V}=3.297636-{V}^{2}$ (29)

Consequently, $V=±\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,n\right)=\left(-29.389788,0.008514\right)$, with normal form coefficient $a=-4.526064$. The eigenvalues are $\left({\lambda }_{1},{\lambda }_{2}\right)=\left(-0.391585,1.96656\right)$. Also, the normal form would be

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

and we have $V=±\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 .

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,n\right)=\left(-11.785736,0.285152\right)$, and parameters value $\left(\varphi ,{I}_{app}\right)=\left(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 :

$\stackrel{˙}{z}=\left(\lambda +i\omega \right)z+{l}_{1}z{|z|}^{2}+{l}_{2}z{|z|}^{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,n\right)=\left(-11.785736,0.285152\right)$ and parameters value $\left(\varphi ,{I}_{app}\right)=\left(0.306345,124.470639\right)$, we have a Bautin point with eigenvalues ${\lambda }_{1,2}=-5.0307+i\left(0.156687\right)$ and the system satisfies the following normal form:

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

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

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

Because, ${l}_{2}<0$, we have supercritical Bautin bifurcation.

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

$\stackrel{˙}{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  .

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,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{˙}{u}=v$

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

where, $a,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,n\right)=\left(-28.744348,0.114090\right)$, $\left(\varphi ,{I}_{app}\right)=\left(0.000000,83.645532\right)$ with $\left(a,b\right)=\left(-5.341083,-1.363867\right)$ and the second Bogdanov-Takens bifurcation happens at $\left(V,n\right)=\left(8.717678,0.610127\right)$, with the parameter values $\left(\varphi ,{I}_{app}\right)=\left(-0.000000,222.452534\right)$ with $\left(a,b\right)=\left(-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 .

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 :

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 :

$\left\{\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+\mathrm{tanh}\left(\left(V-{V}_{1}\right)/{V}_{2}\right)\right],$

${\tau }_{n}\left(V\right)=1/\mathrm{cosh}\left(\left(V-{V}_{3}\right)/\left(2{V}_{4}\right)\right),$

${n}_{\infty }\left(V\right)=\frac{1}{2}\left[1+\mathrm{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 .

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 . 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 .

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: