# Bifurcation Analysis of a Two-Dimensional Neuron Model under Electrical Stimulation

Authors: Chunhua Yuan; Xiangyu Li
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The two-dimensional neuron model can not only reproduce abundant firing patterns, but also satisfy the research of dynamical behavior because of its nonlinear characteristics. It is the most simplified model that includes the fast and slow variables required for neuron firing. In this paper, the dynamic characteristics of two-dimensional neuron model are described by both analytical and numerical methods, and the influence of model parameters and external stimuli on dynamic characteristics is described. The firing characteristics of the Prescott model under external electrical stimulation are studied, and the influence of electrophysiological parameters on the firing characteristics is analyzed. The saddle-node bifurcation and Hopf bifurcation characteristics are studied through the distribution of equilibrium points. It is found that there are critical saddle-node bifurcation and critical Hopf bifurcation in the Prescott model. And the value range of the key parameters that cause the critical bifurcation of the model is obtained by analytical methods.


---

## Body

## 1. Introduction

Neurons are the basic unit of the nervous system, so only by understanding the characteristics and activities of single neurons can we further understand the mystery of neuronal networks and even the operation of the brain [1]. Quantitative mathematical models are an indispensable tool for this purpose [2]. For quantitative models, a large number of detailed descriptions are needed to explain the complex dynamics of a single neuron, and the complexity of the model should be reduced as much as possible while retaining its basic characteristics to achieve the feasibility of model calculations [3]. These two conditions need to achieve a balance.At present, most theoretical research on neurons is based on two types of models, namely, phenomenological models (such as Integrate-and-Fire model [4, 5] and Izhikevich model [6, 7]) and physiological models (such as Hodgkin–Huxley model [8–10] and Morris–Lecar model [11]). The phenomenological models only consider the external input-output relationship and do not consider the internal details. While the detailed physiological models contain abundant biophysical details, it is difficult to describe the dynamic characteristics due to the high dimensions, which is not conducive to calculation and analysis. The two-dimensional differential equations can be studied in an intuitive and visual way through phase plane analysis, which can not only reproduce rich firing patterns, but also satisfy the nonlinear characteristics of dynamic behaviors [12]. It is the most simplified model including the fast and slow variables required for neuronal spiking [13].Neurons are the basic unit of neural information processing, and the generation and conduction of neural information contains rich firing characteristics [14]. If two neurons have the same input, the resulting firing response characteristics are different because of their different dynamic rules. If the membrane conductances of two neurons are different and the input currents are different, they may have the same dynamic rules; therefore they will have the same firing characteristics. Conversely, two neurons with different dynamic rules may have different firing characteristics even if the input currents are the same. Therefore, the characteristics of the neuron dynamic system determine the firing characteristics of the neuron under the given conditions of external input [15–17].Changing the amplitude of the external stimulus will change the phase trajectory and firing state of the neuron. The qualitative change of the phase trajectory of the system is due to the bifurcation process of neuron dynamics [18]. For example, the type of bifurcation determines the excitability of neurons [19]. Chen et al. studied the dynamics of a stochastic Morris–Lecar model of both Type I and II excitability [20]. Li et al. proposed a more realistic hybrid impulsive neuron model based on the Izhikevich model [21]. They analyzed the properties of the equilibrium points and subthreshold bifurcation behavior of the model. Finally, they illustrate the main results of numerical simulations. Bao et al. carried out bifurcation analyses of chaotic and periodic burst firings through constructing the fold and Hopf bifurcation sets of fast spiking subsystem in a 3D autonomous memristor synapse-based Morris–Lecar model [22]. Bao et al. proposed a three-dimensional autonomous Morris–Lecar neuron model to explore periodic bursting behaviors using bifurcation plots, phase plots, and time sequences [23]. Jin et al. proposed a novel class of Caputo-type uncertain random fractional-order model that focuses on the reliability analysis of a fractional-order RC circuit system [24]. However, few researchers studied critical saddle node or critical Hopf bifurcation.Furthermore, the study of the mechanisms of external electrical stimulation on the human brain has become the focus of attention in recent years. On the one hand, external electromagnetic stimulation can change the dynamic behavior of the nervous system. For example, external stimulating currents or external electric fields can affect the excitability of nervous system and the generation and conduction of nervous system information [25–27]. On the other hand, electrical stimulation has become an emerging means of treating neuropsychiatric diseases [28, 29]. More and more bodies of evidence show that regulation and improvement of brain function can be achieved by changing the amplitude, frequency, and other parameters of stimulation [30, 31].The human brain is exposed to an electromagnetic field, which can change the firing characteristics of neurons, which in turn affects neuronal activity. There is an interconversion relationship between stimulation parameters and neuron firing patterns. On the one hand, stimulation parameters can be quantified by different neuron firing patterns. On the other hand, the firing sequence of various patterns can be regarded as the expression of stimulation parameter values. It is the purpose of this paper to confirm that external electrical stimulation can change the dynamic behavior of the nervous system and then electrical stimulation can be used in the field of neurological disease treatment. Therefore, the dynamic characteristics of a two-dimensional neuron model under external electrical stimulation are studied in this paper. The impact of important electrophysiological parameters and external stimulus changes on response firing characteristics, equilibrium point distribution, saddle-node bifurcation, Hopf bifurcation, and critical bifurcation of the neuron model are analyzed by a combination of analytical and numerical methods. The contents of this paper are arranged as follows: a two-dimensional neuron model and its two classes of firing characteristics of excitatory neurons under electrical stimulation are given in Section2. The equilibrium point analysis of the model is given in Section 3. In Sections 4 and 5, the saddle-node bifurcation and Hopf bifurcation at equilibrium points are investigated, respectively. Finally, conclusions and discussions are made in Section 6.


## 2. The Prescott Model and Its Firing Characteristics

The Morris–Lecar model is one of the common biophysical models, which is derived from the experimental study of the electrical activity characteristics of the muscle fibers of the Arctic penguins [32]. Prescott et al. [33] improved the ML model and obtained the dynamic equations composed of a fast variable V and a slow recovery variable w:(1)CmdVdt=I−gLV−EL−gNam∞VV−ENa−gKwV−EK,(2)dwdt=φww∞V−wτwV,where V is the neuron cell membrane voltage, w is the slow ion channel recovery variable, and I is the external stimulating current, gNa, gK, and gL are the maximum conductance of sodium ion channel, the maximum conductance of potassium ion channel, and the leakage conductance, respectively, ENa, EK, and EL are the corresponding reversal potential for sodium, potassium, and chlorine ions, Cm is the cell membrane capacitance of neurons, m∞V is the steady state value of the sodium channel activation variable, w∞V is the steady state value of the potassium channel recovery variable, and τwV is the time constant of the recovery variable. They are all functions of the neuron membrane voltage:(3)m∞V=0.51+tanhV−βmγm,w∞V=0.51+tanhV−βwγw,τwV=1coshV−βw/2γw,where βm and γm are the influencing factors of fast ion channel activation variable and βw and γw are the influencing factors of slow ion channel recovery variable. Changing βw can simulate various firing patterns of neuron models, so βw is selected as the key parameter of the model.The values of the model parameters in this paper are as follows [33]: cm=2μF/cm2, φw=0.15, gL=2mS/cm2, gNa=20mS/cm2, gk=20mS/cm2, EL=−70mV, ENa=50mV, EK=−100mV, βm=−1.2mV, γm=18mV, βw=−10mV, and γw=10mV.Neurons can be divided into two classes according to excitability: neurons of class I excitability can produce firing sequences of any low frequency under different intensities of external stimulation, and theirf−I curves are continuous. Neurons of class II excitability cannot produce firing sequences of any low frequency, and their f−I curves are discontinuous. Repetitive spiking can be produced only when the stimulus intensity reaches above the critical value. The Hodgkin–Huxley model can only reproduce the firing characteristics of class I neurons, and the Morris–Lecar model can only reproduce the firing characteristics of class II neurons. However, in the Prescott model, these two classes of excitability can be reproduced by changing the key parameter βw.Whenβw=−10mV, the firing characteristics are shown in Figure 1 as the input external current increases from I=33.13μA/cm2 to I=34μA/cm2. The neuron model shows no spike, single spike, and periodic spiking, respectively. Neurons can generate arbitrary low frequency firing sequence, so it is Class I excitability at this time. When βw=−20mV, input an external current ranging from I=50μA/cm2 to I=55.9μA/cm2, and its firing characteristics are shown in Figure 2. The neuron model shows no spike at the beginning. Periodic firing only begins when the stimulating current intensity reaches a certain threshold, and the frequency of initial spiking is high, indicating that the neuron cannot produce any low-frequency firing sequence, which represents class II excitability. Generally speaking, class I excitability is caused by saddle-node bifurcation, and class II excitability is caused by Hopf bifurcation. The latter part will analyze these two bifurcation characteristics in detail.Figure 1
Class I firing characteristics of the Prescott model. Whenβw=−10mV, as the input external current increases from I=33.13μA/cm2 to I=34μA/cm2, the neuron model shows no spike, single spike, and periodic spiking, respectively.
(a)(b)(c)(d)Figure 2
Class II firing characteristics of the Prescott model. Whenβw=−20mV, as the input external current increases from I=50μA/cm2 to I=55.9μA/cm2, the neuron model shows no spike at the beginning. Periodic firing only begins when the stimulating current intensity reaches a certain threshold.
(a)(b)(c)(d)


## 3. Equilibrium Point of the Prescott Model

The different firing characteristics of neurons are related to the type and stability of the equilibrium point of the model equation. The equilibrium point of the Prescott model can be determined by its dynamic equations. LetdV/dt=0 and dw/dt=0 in formulas (1) and (2); we can get the corresponding zero line equations:(4)w=−gLV−EL−0.5gNa1+tanhV−βm/γmV−ENa+IgKV−EK,(5)w=w∞V=0.51+tanhV−βwγw.When the two zero lines intersect, that is, whenw satisfies formulas (4) and (5) at the same time, the equilibrium point of the system equation is obtained. The corresponding membrane voltage V and recovery variable w are the equilibrium point V0,w0 of the Prescott model.Let(6)f=I−gLV−EL−gNam∞VV−ENa−gKwV−EK.Then,f is the total current of neuron cell membrane:(7)f=I−gLV−EL−0.5gNa1+tanhV−βmγmV−ENa−0.5gK1+tanhV−βwγwV−EK.Let the external currentI be the variable parameter to draw the f−V curve (V is the abscissa, f is the ordinate). The intersection point of the curve with the horizontal axis is the equilibrium point V0,w0 of the neuron model.The change of the electrophysiological parameters of neurons leads to the change of the shape and relative position of the zero line of the model equation, thus changing the number and type of the equilibrium point of the equation. In the Prescott model, the key parameterβw mainly affects the position of the zero line w, as shown in Figure 3. When βw decreases, w zero line shifts to the left, and when βw increases, w zero line shifts to the right. The left and right translation of the w zero line changes the position of its intersection with the V zero line, thereby changing the membrane voltage at the equilibrium point. With the change of βw=−20mV∼0mV, the type of equilibrium point changes from a single focus to the coexistence of a focus, a saddle, and a node, and finally evolves into a single node.Figure 3
The evolution process of equilibrium points in the Prescott model at differentβw: (a) total current of neuron cell membrane; (b) phase plane.
(a)(b)


## 4. Analysis of Saddle-Node Bifurcation at Equilibrium Point

With the change of the bifurcation parameters, the stable node and the unstable saddle gradually approach. When the bifurcation parameter is equal to the bifurcation value, the saddle and the node collide with each other and disappear, and the saddle-node bifurcation occurs. The newly generated equilibrium point is neither a saddle nor a node. It is called a saddle node, which is characterized by the fact that one of the two real eigenvalues is zero. The saddle node is stable in half of its neighborhood, and unstable in the other half of its neighborhood [18].


### 4.1. Determination of the Saddle Node

For neurons, a stable equilibrium point corresponds to the resting state. When the stable equilibrium point and the unstable equilibrium point collide and disappear, the resting state will disappear, and the phase trajectory of the system will become a stable limit cycle. The neuron begins to fire periodically. As the stimulus intensity increases, theV zero line moves upward. When the stimulus intensity reaches the threshold (the minimum stimulus intensity that can cause a spike), the stable node among the original three intersection points collides with the unstable saddle and disappears. A stable limit cycle is generated and the neuron produces repetitive firing.It can be seen from the above that when theV zero line of the neuron model is tangent to the w zero line, that is, when the total current f is tangent to the horizontal axis, saddle-node bifurcation will occur. Therefore, df/dV=0 at the saddle node:(8)dfdV=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞′VV−EK−gKw∞V=0.After obtaining the membrane potentialVSN at the saddle node, letting f=0 in formula (7) can determine the external stimulating current ISN when the saddle-node bifurcation occurs:(9)I=gLV−EL+0.5gNa1+tanhV−βmγmV−ENa+0.5gK1+tanhV−βwγwV−EK.Figure4 shows the derivative curve df/dV of the right function of the model equation. The membrane potential at the saddle node can be obtained from f˙V=0: VSN1=−38.8483mV, VSN2=−33.7954mV. From formula (9), the stimulating current at the saddle node can be obtained as follows: ISN1=33.1855μA/cm2, ISN2=32.7822μA/cm2.Figure 4
Thedf/dV curve of the Prescott model.WhenISN1=33.1855μA/cm2, the total current f and the phase plane of the neuron model are shown in Figure 5. f is tangent to the horizontal axis at VSN1=−38.8483mV, where the saddle and the node collide and disappear, and a saddle-node bifurcation occurs. When ISN2=32.7822μA/cm2, the total current f and the phase plane of the neuron model are shown in Figure 6. f is tangent to the horizontal axis at VSN2=−33.7954mV, where the saddle and the node collide and disappear, and another saddle-node bifurcation occurs.Figure 5
Saddle-node bifurcation occurs in the Prescott model whenI=33.1855μA/cm2: (a) total current of cell membrane; (b) phase plane.
(a)(b)Figure 6
Saddle-node bifurcation occurs in the Prescott model whenI=32.7822μA/cm2: (a) total current of cell membrane; (b) phase plane.
(a)(b)The characteristic coefficient distributionq−p corresponding to the equilibrium point of the Prescott model is shown in Figure 7, where p is the first order coefficient of the characteristic equation of the equilibrium point and q is the constant term of the characteristic equation of the equilibrium point. In the figure, the black “+” line, the blue “×” line, and the pink “o” line represent the node, the saddle, and the focus, respectively. The blue diamond represents the saddle node, and the red curve represents p2−4q=0, which is the dividing line between focus and node of the equilibrium point. It can be seen from the figure that the blue curve q<0, so it is the saddle. The black curve satisfies p2−4q>0, so it is the node. The pink curve satisfies p2−4q<0, so it is the focus. There is a very small distance of node between the pink curve and the blue curve. Therefore, the neuron model has two saddle-node bifurcations at ISN1=33.1855μA/cm2 and ISN2=32.7822μA/cm2.Figure 7
Distribution of characteristic coefficientsq−p of equilibrium points of the Prescott model. The red line represents p2−4q=0 curve.


### 4.2. Multivalued Equilibrium Point Region and Critical Saddle Node

From the analysis and calculation in the previous section, it can be seen that the derivative of the total current is not a monotonic function. There are two tangent points between the total currentf and the horizontal axis, indicating that the neuron model has two saddle nodes. Let ΔVSN=VSN2−VSN1 be the difference between the membrane potentials of the two saddle nodes and ΔISN=ISN1−ISN2 be the difference between the external currents corresponding to the two saddle nodes, as shown in Figure 8.Figure 8
The difference between the membrane potentials and the corresponding external currents of the two saddle nodes.When the external current is betweenISN1 and ISN2, the number of equilibrium points of the neuron model is not unique, which means it is a multivalued equilibrium point region. When βw changes, the range of the multivalued equilibrium point region also changes. As βw decreases, ΔVSN and ΔISN gradually decrease until they equal zero, and the multivalued equilibrium point region disappears. It can be seen from Figure 9 that ΔV and ΔI show a nonlinear decrease as βw decreases. When βw=−10.4176mV, ΔVSN⟶0, and ΔISN⟶0, the two saddle nodes meet and disappear. At this time, regardless of the value of the external current, the neuron model has only one equilibrium point, and no more multivalued equilibrium point region is generated. The boundary point of the single and multivalued equilibrium point region is called the critical saddle node in this paper.Figure 9
The change of (a) the difference of membrane potentials and (b) the difference of external currents withβw of two saddle nodes.
(a)(b)The total current and phase plane of the Prescott model at the critical saddle node are shown in Figure10. At this time, the membrane potential is VSN0=36.66mV, and the corresponding external stimulating current is ISN0=33.57μA/cm2.Figure 10
The (a) total current and (b) phase plane at the critical saddle node.
(a)(b)The value ofβw at the critical saddle node is calculated below. At the critical saddle node, the total current f is tangent to the horizontal axis and there is only one tangent point, as shown in Figure 10(a), so the second derivative of f is equal to zero; that is, df/dV=d2f/dV2=0:(10)d2fdV2=−gNadm∞′VdVV−ENa−2gNam∞′V−gKdw∞′VdVV−EK−2gKw∞′V=0.Therefore,(11)−gNadm∞′VdVV−ENa−2gNam∞′V−gKdw∞′VdVV−EK−2gKw∞′V=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞′VV−EK−gKw∞V.After calculation, we get(12)−gNashV−βm/γmγm2ch3V−βm/γmV−ENa+0.5γm1ch2V−βm/γmV−ENa−2+m∞V=gKshV−βm/γmγw2ch3V−βm/γmV−EK+0.5γw1ch2V−βm/γmV−EK−2+w∞V+gL.Thus, the calculation formula ofβw corresponding to the critical saddle node can be obtained. Figure 11(a) shows the curves of the first derivative and the second derivative of f when βw=−10.4176mV, which correspond to the critical saddle node. Figure 11(b) shows the curves of the first derivative and the second derivative of f when βw=−15mV. It can be seen from the figure that df/dV=d2f/dV2=0 when the neuron model is at the critical saddle node. While βw<−10.4176mV, d2f/dV2≠0 when df/dV=0, the neuron model will no longer have saddle-node bifurcation and there will be no multivalued equilibrium point region.Figure 11
The first and the second derivative of the total currentf when βw is different: (a)βw=−10.4176mV; (b)βw=−15mV.
(a)(b)


## 4.1. Determination of the Saddle Node

For neurons, a stable equilibrium point corresponds to the resting state. When the stable equilibrium point and the unstable equilibrium point collide and disappear, the resting state will disappear, and the phase trajectory of the system will become a stable limit cycle. The neuron begins to fire periodically. As the stimulus intensity increases, theV zero line moves upward. When the stimulus intensity reaches the threshold (the minimum stimulus intensity that can cause a spike), the stable node among the original three intersection points collides with the unstable saddle and disappears. A stable limit cycle is generated and the neuron produces repetitive firing.It can be seen from the above that when theV zero line of the neuron model is tangent to the w zero line, that is, when the total current f is tangent to the horizontal axis, saddle-node bifurcation will occur. Therefore, df/dV=0 at the saddle node:(8)dfdV=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞′VV−EK−gKw∞V=0.After obtaining the membrane potentialVSN at the saddle node, letting f=0 in formula (7) can determine the external stimulating current ISN when the saddle-node bifurcation occurs:(9)I=gLV−EL+0.5gNa1+tanhV−βmγmV−ENa+0.5gK1+tanhV−βwγwV−EK.Figure4 shows the derivative curve df/dV of the right function of the model equation. The membrane potential at the saddle node can be obtained from f˙V=0: VSN1=−38.8483mV, VSN2=−33.7954mV. From formula (9), the stimulating current at the saddle node can be obtained as follows: ISN1=33.1855μA/cm2, ISN2=32.7822μA/cm2.Figure 4
Thedf/dV curve of the Prescott model.WhenISN1=33.1855μA/cm2, the total current f and the phase plane of the neuron model are shown in Figure 5. f is tangent to the horizontal axis at VSN1=−38.8483mV, where the saddle and the node collide and disappear, and a saddle-node bifurcation occurs. When ISN2=32.7822μA/cm2, the total current f and the phase plane of the neuron model are shown in Figure 6. f is tangent to the horizontal axis at VSN2=−33.7954mV, where the saddle and the node collide and disappear, and another saddle-node bifurcation occurs.Figure 5
Saddle-node bifurcation occurs in the Prescott model whenI=33.1855μA/cm2: (a) total current of cell membrane; (b) phase plane.
(a)(b)Figure 6
Saddle-node bifurcation occurs in the Prescott model whenI=32.7822μA/cm2: (a) total current of cell membrane; (b) phase plane.
(a)(b)The characteristic coefficient distributionq−p corresponding to the equilibrium point of the Prescott model is shown in Figure 7, where p is the first order coefficient of the characteristic equation of the equilibrium point and q is the constant term of the characteristic equation of the equilibrium point. In the figure, the black “+” line, the blue “×” line, and the pink “o” line represent the node, the saddle, and the focus, respectively. The blue diamond represents the saddle node, and the red curve represents p2−4q=0, which is the dividing line between focus and node of the equilibrium point. It can be seen from the figure that the blue curve q<0, so it is the saddle. The black curve satisfies p2−4q>0, so it is the node. The pink curve satisfies p2−4q<0, so it is the focus. There is a very small distance of node between the pink curve and the blue curve. Therefore, the neuron model has two saddle-node bifurcations at ISN1=33.1855μA/cm2 and ISN2=32.7822μA/cm2.Figure 7
Distribution of characteristic coefficientsq−p of equilibrium points of the Prescott model. The red line represents p2−4q=0 curve.


## 4.2. Multivalued Equilibrium Point Region and Critical Saddle Node

From the analysis and calculation in the previous section, it can be seen that the derivative of the total current is not a monotonic function. There are two tangent points between the total currentf and the horizontal axis, indicating that the neuron model has two saddle nodes. Let ΔVSN=VSN2−VSN1 be the difference between the membrane potentials of the two saddle nodes and ΔISN=ISN1−ISN2 be the difference between the external currents corresponding to the two saddle nodes, as shown in Figure 8.Figure 8
The difference between the membrane potentials and the corresponding external currents of the two saddle nodes.When the external current is betweenISN1 and ISN2, the number of equilibrium points of the neuron model is not unique, which means it is a multivalued equilibrium point region. When βw changes, the range of the multivalued equilibrium point region also changes. As βw decreases, ΔVSN and ΔISN gradually decrease until they equal zero, and the multivalued equilibrium point region disappears. It can be seen from Figure 9 that ΔV and ΔI show a nonlinear decrease as βw decreases. When βw=−10.4176mV, ΔVSN⟶0, and ΔISN⟶0, the two saddle nodes meet and disappear. At this time, regardless of the value of the external current, the neuron model has only one equilibrium point, and no more multivalued equilibrium point region is generated. The boundary point of the single and multivalued equilibrium point region is called the critical saddle node in this paper.Figure 9
The change of (a) the difference of membrane potentials and (b) the difference of external currents withβw of two saddle nodes.
(a)(b)The total current and phase plane of the Prescott model at the critical saddle node are shown in Figure10. At this time, the membrane potential is VSN0=36.66mV, and the corresponding external stimulating current is ISN0=33.57μA/cm2.Figure 10
The (a) total current and (b) phase plane at the critical saddle node.
(a)(b)The value ofβw at the critical saddle node is calculated below. At the critical saddle node, the total current f is tangent to the horizontal axis and there is only one tangent point, as shown in Figure 10(a), so the second derivative of f is equal to zero; that is, df/dV=d2f/dV2=0:(10)d2fdV2=−gNadm∞′VdVV−ENa−2gNam∞′V−gKdw∞′VdVV−EK−2gKw∞′V=0.Therefore,(11)−gNadm∞′VdVV−ENa−2gNam∞′V−gKdw∞′VdVV−EK−2gKw∞′V=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞′VV−EK−gKw∞V.After calculation, we get(12)−gNashV−βm/γmγm2ch3V−βm/γmV−ENa+0.5γm1ch2V−βm/γmV−ENa−2+m∞V=gKshV−βm/γmγw2ch3V−βm/γmV−EK+0.5γw1ch2V−βm/γmV−EK−2+w∞V+gL.Thus, the calculation formula ofβw corresponding to the critical saddle node can be obtained. Figure 11(a) shows the curves of the first derivative and the second derivative of f when βw=−10.4176mV, which correspond to the critical saddle node. Figure 11(b) shows the curves of the first derivative and the second derivative of f when βw=−15mV. It can be seen from the figure that df/dV=d2f/dV2=0 when the neuron model is at the critical saddle node. While βw<−10.4176mV, d2f/dV2≠0 when df/dV=0, the neuron model will no longer have saddle-node bifurcation and there will be no multivalued equilibrium point region.Figure 11
The first and the second derivative of the total currentf when βw is different: (a)βw=−10.4176mV; (b)βw=−15mV.
(a)(b)


## 5. Analysis of Hopf Bifurcation at Equilibrium Point

As the bifurcation parameters change, the stability of the equilibrium point of the nonlinear system changes. The stable focus turns into an unstable focus and a limit cycle is generated near it, and the Hopf bifurcation phenomenon occurs. At this time, the equilibrium point becomes the center point, and its eigenvalues become pure imaginary number. When the bifurcation parameter changes to the bifurcation value, the unstable limit cycle shrinks to a stable equilibrium point and turns it into an unstable equilibrium point.


### 5.1. Hopf Bifurcation of Two-Dimensional Prescott Model

For a two-dimensional nonlinear system, the necessary prerequisite for Hopf bifurcation is that the corresponding linear system has center points. The two eigenvalues of the equilibrium point of a two-dimensional system are(13)λ1,2=12−p±p2−4q.Whenp=0 and q>0, the eigenvalues λ1,2=±−q are pure imaginary numbers, and the equilibrium point is center point. Therefore, the conditions for the Hopf bifurcation of the system are obtained: (1) p=0; and (2) q>0.For the Prescott model:(1)
p=0, in the Prescott model:(14)p=gL+gNam∞′VV−ENa+gNam∞V+gKwCm+φwcoshV−βw2γw=0.Therefore,(15)w=−CmφwcoshV−βw/2γw−gL−gNam∞′VV−ENa−gNam∞VgK.If the equilibrium pointV0,W0 of the model satisfies formula (15), then the coefficient of the first-order term of the characteristic equation of the equilibrium point is zero.(2)
q>0:According to condition (1),p=a+d=0, so q=ad−bc=−a2−bc>0, that is, a2+bc<0. Therefore, this condition is equivalent to ① bc<0, that is, b and c are opposite signs, and ② bc>a2 or bc>d2.In the Prescott model,b=−gKV−EK/Cm and EK=−100mV. Usually, the membrane voltage is not less than −100mV, so V>EK. To command V>EK at this time, only let c > 0; that is,(16)c=φww∞′VτwV−w∞V−wτw′VτwV2>0.When the system is at equilibrium,w=w∞V, thus(17)φww∞′VτwV>0.Substitutingw∞′V=0.5/γw1/cosh2V−βw/γw, τwV=1/coshV−βw/2γw, we have(18)φw2γw1cosh2V−βw/γw⋅coshV−βw/2γwis always greater than0.That is, when the neuron model equation is under given parameters and the coefficient of its characteristic equationp=0, condition ① is automatically satisfied.Condition ②:bc>d2, that is,(19)−gKV−EKφww∞′VCmτwV>φw2cosh2V−βw2γw.The hyperbolic cosine function is a constant positive function, so we can get(20)−gKV−EK2Cmγwφwcosh2V−βw/γwcoshV−βw/γw>1.Therefore, when the equilibrium point of the Prescott model satisfies both formulas (15) and (20), Hopf bifurcation occurs in the system.


### 5.2. Critical Hopf Bifurcation

The key parameterβw has an important influence on Hopf bifurcation. When βw is different, the characteristic equation coefficients of the equilibrium point are also different. Only when p=0 and q>0 can the model produce Hopf bifurcation.Whenβw is different, the characteristic equation coefficient q−p of the equilibrium point is shown in Figure 12. It can be seen from the figure that when βw=−25mV, there is no intersection point between the p curve and the vertical axis. At this time, the model will not have Hopf bifurcation. And when βw<−25mV, the value of p will no longer be equal to zero, so the model will not have Hopf bifurcation when βw<−25mV. When βw=0mV and βw=−5mV, if p=0, q<0, the model will not have Hopf bifurcation. Therefore, Hopf bifurcation occurs only when the value of βw is within a certain range. The boundary of which the model will have no Hopf bifurcation is called critical Hopf bifurcation in this article. The value of βw at the critical Hopf bifurcation point is calculated below.Figure 12
Theq−p curve of the characteristic equation of the equilibrium point when βw is different. When βw≤−25mV, there is no intersection point between the p curve and the vertical axis. At this time, the model will not have Hopf bifurcation. When βw≥−5mV, if p=0, then q<0, the model will not have Hopf bifurcation.Through the above qualitative analysis, we can see that there are two critical Hopf bifurcation points. Whenp=0, q<0, the model no longer has Hopf bifurcation. In this case, this boundary point is called the upper critical Hopf bifurcation point. Therefore, the critical condition is when p=0, q=0; that is,(21)p=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞Vcm−φwcoshV−βw2γw=0,q=−φw2cosh2V−βw2γw+gKV−EKφww∞′VcmτwV=0.The value ofβw and the membrane voltage at the upper critical Hopf bifurcation point can be obtained from above. Figure 13 depicts the variation curves of characteristic coefficients p and q with membrane voltage V. The pink line is the p−V curve, and the blue line is the q−V curve. When βw1=−9.8461mV and VHo1=−39.1726mV, the model satisfies p=0 and q=0; hence it is the upper critical Hopf bifurcation point of the model.Figure 13
p,q−V curves at the upper critical Hopf bifurcation point.When thep curve no longer intersects the zero axis, the model no longer has Hopf bifurcation. This boundary point is called the lower critical Hopf bifurcation point. Therefore, the critical condition is that the p−V curve is tangent to the zero axis; that is,(22)p=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞VCm−φwcoshV−βw2γw=0,dpdV=−−gNam∞″VV−ENa−2gNam∞′V−gKw∞′VCm−φw2γwshV−βw2γw=0,where m∞″V=−shV−βm/γm/γm2ch3V−βm/γm, m∞′V=0.5/γm1/ch2V−βm/γm, w∞′V=0.5/γw1/ch2V−βw/γw.The value ofβw and membrane voltage at the lower critical Hopf bifurcation point are calculated as βw2=−22.9899mV, VHo2=−32.52mV, and the corresponding input current is IHo=180.0713μA/cm2. The p curve, dp/dV curve, q curve, and input current I curve at the lower critical bifurcation point are shown in Figure 14, and the I curve in the figure is drawn at a ratio of 1%.Figure 14
p,q−V curves at the lower critical Hopf bifurcation point.Whenβw1=−9.8461mV and βw2=−22.9899mV, the characteristic equation coefficient q−p curves corresponding to the upper and lower critical Hopf bifurcations of the model are shown in Figure 15. The value of βw corresponding to the area between the two curves is the range of the value of βw where Hopf bifurcation can occur; that is, −22.9899mV<βw<−9.8461mV.Figure 15
Characteristic coefficientq−p curves corresponding to two critical Hopf bifurcations.


## 5.1. Hopf Bifurcation of Two-Dimensional Prescott Model

For a two-dimensional nonlinear system, the necessary prerequisite for Hopf bifurcation is that the corresponding linear system has center points. The two eigenvalues of the equilibrium point of a two-dimensional system are(13)λ1,2=12−p±p2−4q.Whenp=0 and q>0, the eigenvalues λ1,2=±−q are pure imaginary numbers, and the equilibrium point is center point. Therefore, the conditions for the Hopf bifurcation of the system are obtained: (1) p=0; and (2) q>0.For the Prescott model:(1)
p=0, in the Prescott model:(14)p=gL+gNam∞′VV−ENa+gNam∞V+gKwCm+φwcoshV−βw2γw=0.Therefore,(15)w=−CmφwcoshV−βw/2γw−gL−gNam∞′VV−ENa−gNam∞VgK.If the equilibrium pointV0,W0 of the model satisfies formula (15), then the coefficient of the first-order term of the characteristic equation of the equilibrium point is zero.(2)
q>0:According to condition (1),p=a+d=0, so q=ad−bc=−a2−bc>0, that is, a2+bc<0. Therefore, this condition is equivalent to ① bc<0, that is, b and c are opposite signs, and ② bc>a2 or bc>d2.In the Prescott model,b=−gKV−EK/Cm and EK=−100mV. Usually, the membrane voltage is not less than −100mV, so V>EK. To command V>EK at this time, only let c > 0; that is,(16)c=φww∞′VτwV−w∞V−wτw′VτwV2>0.When the system is at equilibrium,w=w∞V, thus(17)φww∞′VτwV>0.Substitutingw∞′V=0.5/γw1/cosh2V−βw/γw, τwV=1/coshV−βw/2γw, we have(18)φw2γw1cosh2V−βw/γw⋅coshV−βw/2γwis always greater than0.That is, when the neuron model equation is under given parameters and the coefficient of its characteristic equationp=0, condition ① is automatically satisfied.Condition ②:bc>d2, that is,(19)−gKV−EKφww∞′VCmτwV>φw2cosh2V−βw2γw.The hyperbolic cosine function is a constant positive function, so we can get(20)−gKV−EK2Cmγwφwcosh2V−βw/γwcoshV−βw/γw>1.Therefore, when the equilibrium point of the Prescott model satisfies both formulas (15) and (20), Hopf bifurcation occurs in the system.


## 5.2. Critical Hopf Bifurcation

The key parameterβw has an important influence on Hopf bifurcation. When βw is different, the characteristic equation coefficients of the equilibrium point are also different. Only when p=0 and q>0 can the model produce Hopf bifurcation.Whenβw is different, the characteristic equation coefficient q−p of the equilibrium point is shown in Figure 12. It can be seen from the figure that when βw=−25mV, there is no intersection point between the p curve and the vertical axis. At this time, the model will not have Hopf bifurcation. And when βw<−25mV, the value of p will no longer be equal to zero, so the model will not have Hopf bifurcation when βw<−25mV. When βw=0mV and βw=−5mV, if p=0, q<0, the model will not have Hopf bifurcation. Therefore, Hopf bifurcation occurs only when the value of βw is within a certain range. The boundary of which the model will have no Hopf bifurcation is called critical Hopf bifurcation in this article. The value of βw at the critical Hopf bifurcation point is calculated below.Figure 12
Theq−p curve of the characteristic equation of the equilibrium point when βw is different. When βw≤−25mV, there is no intersection point between the p curve and the vertical axis. At this time, the model will not have Hopf bifurcation. When βw≥−5mV, if p=0, then q<0, the model will not have Hopf bifurcation.Through the above qualitative analysis, we can see that there are two critical Hopf bifurcation points. Whenp=0, q<0, the model no longer has Hopf bifurcation. In this case, this boundary point is called the upper critical Hopf bifurcation point. Therefore, the critical condition is when p=0, q=0; that is,(21)p=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞Vcm−φwcoshV−βw2γw=0,q=−φw2cosh2V−βw2γw+gKV−EKφww∞′VcmτwV=0.The value ofβw and the membrane voltage at the upper critical Hopf bifurcation point can be obtained from above. Figure 13 depicts the variation curves of characteristic coefficients p and q with membrane voltage V. The pink line is the p−V curve, and the blue line is the q−V curve. When βw1=−9.8461mV and VHo1=−39.1726mV, the model satisfies p=0 and q=0; hence it is the upper critical Hopf bifurcation point of the model.Figure 13
p,q−V curves at the upper critical Hopf bifurcation point.When thep curve no longer intersects the zero axis, the model no longer has Hopf bifurcation. This boundary point is called the lower critical Hopf bifurcation point. Therefore, the critical condition is that the p−V curve is tangent to the zero axis; that is,(22)p=−gL−gNam∞′VV−ENa−gNam∞V−gKw∞VCm−φwcoshV−βw2γw=0,dpdV=−−gNam∞″VV−ENa−2gNam∞′V−gKw∞′VCm−φw2γwshV−βw2γw=0,where m∞″V=−shV−βm/γm/γm2ch3V−βm/γm, m∞′V=0.5/γm1/ch2V−βm/γm, w∞′V=0.5/γw1/ch2V−βw/γw.The value ofβw and membrane voltage at the lower critical Hopf bifurcation point are calculated as βw2=−22.9899mV, VHo2=−32.52mV, and the corresponding input current is IHo=180.0713μA/cm2. The p curve, dp/dV curve, q curve, and input current I curve at the lower critical bifurcation point are shown in Figure 14, and the I curve in the figure is drawn at a ratio of 1%.Figure 14
p,q−V curves at the lower critical Hopf bifurcation point.Whenβw1=−9.8461mV and βw2=−22.9899mV, the characteristic equation coefficient q−p curves corresponding to the upper and lower critical Hopf bifurcations of the model are shown in Figure 15. The value of βw corresponding to the area between the two curves is the range of the value of βw where Hopf bifurcation can occur; that is, −22.9899mV<βw<−9.8461mV.Figure 15
Characteristic coefficientq−p curves corresponding to two critical Hopf bifurcations.


## 6. Conclusions

This paper uses both analytical and numerical methods to study the dynamic characteristics of a two-dimensional neuron model and uses the changes of key parameters to describe the relationship between the input and output characteristics of the neuron and the internal parameters and input stimuli. The firing characteristics of the two-dimensional Prescott model under external electrical stimulation are studied. Through the analysis of the equilibrium point of the model, it is found that the electrophysiological parameters have an important influence on the distribution and bifurcation of the equilibrium point. The critical saddle node of the Prescott model is related to the key parametersβw. When βw<−10.4176mV, the saddle-node bifurcation disappears. At this time, regardless of the value of the external current, the model has only one equilibrium point and no more multivalued equilibrium point region is generated. The Prescott model also has critical Hopf bifurcation. Only when βw is within a certain parameter range will the model have Hopf bifurcation, and there are upper and lower critical values. When −22.9899mV<βw<−9.8461mV, the model can undergo Hopf bifurcation.In saddle-node bifurcation, when external electrical stimuli of different intensities are input, the saddle and the node approach each other, then collide, and disappear, resulting in a stable limit cycle, which produces a periodic firing that has a continuousf−I curve. In the Hopf bifurcation, the stable equilibrium point loses its stability and gives birth to a stable limit cycle, resulting in a periodic firing mode of the discontinuous f−I curve. This paper analyzes the saddle-node bifurcation and Hopf bifurcation of the Prescott model under electrical stimulation. In particular, the critical saddle-node bifurcation and the critical Hopf bifurcation are obtained by the combination of numerical analysis and analytical methods. These two bifurcations correspond to different firing patterns of neurons, and the microscopic basis of many neuropsychiatric diseases is the abnormal firing of brain neurons, which in turn leads to changes in neural information coding. Therefore, the results of this article are helpful to the research of electromagnetic field in the neurological disease treatment to change the firing rhythms of morbid neurons. And it is also of great significance in the research of nerve control and nerve electrical stimulation, which can provide a certain theoretical basis for understanding the influence of electromagnetic fields on brain neurons.



---

*Source: 9314736-2021-09-17.xml*