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

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 de-scribe the phenomenon of neural bursting, and investigate the dynamics of Morris-Lecar model as a square-wave burster, elliptic burster and parabolic burster.


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.

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 where V demonstrates membrane potential, and n the activation variable of the persistent K + current, so it is a two-dimensional vector ( ) 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 Ca K − are the major ionic currents. Therefore where Ca I (inward current) is negative and also K I 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, app I has not been included since we consider it as a bifurcation parameter for bifurcation diagrams.

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  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 4 displays the occurence of limit cycle corresponding to hopf bifurcation with increasing app I .
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 20 mv V = − . Also, depends on the initial values of membrane potential, the size of the excursion is different.
Using continuation software Matcont and choosing app I as a bifurcation parameter, we could detect numerically hopf bifurcation points. For    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 ( ) ( ) . Since, the real part is negative, it means that the equilibrium point is stable. The topological normal form for hopf bifurcation has the form: Here, β does not have any dynamical effect. The normal form for first hopf has the form below: Here, 0 θ >  means that the frequency of damped or sustained oscillations around this point, 0 ω , is positive and increasing. But for the othe normal form:

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 app I .    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 3 6.13675 5.317042 . The normal form for this bifurcation can be written as: and for this case, we can write the following normal form:

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)   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 1 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   Because here 1 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.

Co-Dimension Two Bifurcations
In this section, we focus on co-dimension two bifurcations with app I and φ 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 2 l are non zero, Because, 2 0 l < , we have supercritical Bautin bifurcation.
Generalized hopf (Bautin) bifurcation in polar coordinates has the following normal form [14]: In our simulation, the curve LPC corresponds to the saddle-node bifurcation of periodic orbits. As we can see in Figure 10, for 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: Because of these two conditions, Bogdanov-Takens is a codimension two bifurcation that has the following normal form u v =  . 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.

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

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 app I 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]:  The parameter Ca K implies to the calcium removal rate and ε represents the ratio of free to total calcium in the cell. Because calcium is highly buffered, ε 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, KCa I called outward current. If conductance KCa g z is large, the cell has hyper-polarization state which is corresponding to resting behavior. Conversely, if KCa g 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 app I in Figure 15.

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 app I in Figure 16.

Morris-Lecar Model as a Parabolic Burster
Unlike two previous bursters which we need only one slow variable for bursting   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 be-cause of the bistability and hysteresis loop. In parabolic burster, the model has the form: Ca 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 app I in Figure 18.

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 app I and temperature φ . 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.