Firing patterns and bifurcation analysis of neurons under electromagnetic induction

1.   Introduction

The electrical activities of neurons are the basis of transmitting information in the nervous system. Especially, the excitability of neurons is sensitive to many factors, such as the neuronal intracellular and extracellular ion concentration, environmental noise, temperature and so on. The single neuron can exhibit multiple firing patterns under the external stimulus, while the transformation between different firing patterns corresponds to the bifurcation of models, such as the pancreas $\beta$ cells increase the duty cycle as blood sugar levels rise [17]. Duan et al. study the bifurcation associated with the transition of bursting patterns in the Chay model [5]. MMOs corresponding to more complex bifurcations have also been extensively investigated through geometric singular perturbation theory [25]. Teka Krasimira et al. discuss in detail the transition from plateau to pseudo-plateau bursting by using two-parameter bifurcation [22]. It can be seen that the bifurcation theory is a powerful tool to analyze the firing patterns of neurons.

According to Faraday's law of electromagnetic induction, the action potential of excitable neurons can generate the magnetic field in the medium, and the magnetic field can influence the neurons in turn. Interestingly, Barry et al. have found the magnetic field, indicating that the magnetic field exists generated by action potential [1]. The exogenous magnetic field has obvious influence on the electrical activity of neurons [20][21]. Therefore, it is a direction to study theoretically the dynamic behavior of neurons under electromagnetic induction. Recently, many researchers have begun to put attention to the dynamics of electromagnetic induction on neurons and neural networks [4]. Li et al. build a mathematical model under electromagnetic radiation to study the variation of firing rhythms of neurons [11]. Then, Lv and Ma introduce a magnetic flux variable and the feed-back current on the membrane potential on the Hindmarsh-Rose (HR) neuron model and study the electrical activities [14]. On this basis, Mvogo et al. investigate the spatiotemporal formation of patterns in a diffusive FitzHugh-Nagumo network under the effects of electromagnetic induction [16]. Further, most investigations focus on the simplified neuronal models [15][9], while the effects of electromagnetic induction on the biological neuron models are relatively little considered.

Bursting is characteristic of the nervous system. Bursting oscillations play an important role in information transmission of the nervous system. For this reason, many models have been established to study different bursting activities [12]. Studies show that bursting consists of a fast process and a slow process: the generation of action potentials and modulates fast variables [18]. Thus, the fast process is separated from the slow process, then we can analyze the firing behavior of the fast subsystem, which is called fast-slow dynamics analysis [18]. Izhikevich uses this method to topologically classify bursting types [8]. Moreover, Bertram has discussed the application of fast-slow dynamics analysis in different oscillations in detail, such as relaxation oscillations, neuronal bursting oscillations, canard oscillations, and MMOs [2].

In this paper, taking the endocrine model [23] as the research object, we aim to explore the bifurcation and firing rhythms of neurons under electromagnetic induction and external forcing direct current. The organization of this paper is as follows. In Section 2, we describe the endocrine model and introduce electromagnetic induction. In Section 3, for one-parameter bifurcation, we study the codimension-one and Interspike Intervals bifurcation, and the first Lyapunov coefficient is computed to judge the Hopf bifurcation stability. At the same time, the chaotic region of the system is discussed through the first and second Lyapunov exponents. In Section 4, according to the difference of coupling strength between membrane potential of the neuron and magnetic flux, the firing behaviors of the neuron model are discussed. In Section 5, we also explore two-parameter bifurcation analysis in the $( I_{ext},k_0\ )$ phase plane, and investigate the Bogdanov-Takens bifurcation. Furthermore, we get three bifurcation curves: a saddle-node bifurcation curve, a Hopf bifurcation curve and a homoclinic bifurcation curve.

2.   Model description

To study the plateau-bursting patterns, Tsaneva-Atanasova et al. have constructed a generic simplified endocrine model based on elements drawn from several published models [23]. Here we introduce magnetic flux across the membrane and feedback current as induction current resulted from the variation of magnetic flux and field, then we get the following system:

where variable $V$ is the membrane potential, $n$ is the gating variables for the voltage-gated potassium, $c$ represents the free ${Ca}^{2+}$ concentration in the cytosol the corresponding fluxes and $\varphi$ is magnetic flux across the membrane. The ionic currents $I_{Ca}$, $I_K$ and $I_{K_{Ca}}$ represent inward calcium current, delated-rectifier potassium current and calcium-activated potassium ion current, respectively. $I_{ext}$ is the external forcing direct current. $C_m = {10}^{-5}\times A_{cell}$ is the membrane capacitance; $\tau_n$ is the activation time constant for the delayed rectifier channel; $n_\infty$ is the steady-state function for gating variable $n$; $f_c$ is the fraction of free to total cytosolic ${Ca}^{2+}$; $\theta = {10}^5\left(2\times9.65\times A_{cell}\right)^{-1}$ is a factor that converts current to flux, where $A_{cell} = \pi\times d_{cell}^2$ is the surface area of cell; $k_{PMCA}$ is the plasma membrane ${Ca}^{2+}$ ATPase pump rate. The expressions of ion current are as following:

The steady-state activation functions are described as the following equations:

In addition, other fixed parameters used in this paper are given in Table 1.

In this paper, the term of the $k_0V\rho\left(\varphi\right)$ describes the feedback current as induction current results from the variation of magnetic flux and $k_0$ is feedback gain. $\rho\left(\varphi\right) = \alpha+3\beta\varphi^2$ is a magnetic flux-controlled memristor, which is equivalent to memory conductance and is used to describe the coupling relation between membrane potential of the neuron and magnetic flux, $\alpha$ and $\beta$ are fixed parameter values [14]. Parameter $k_1$ and $k_2$ are the interaction between membrane potential and magnetic flux [14].

3.   One-parameter bifurcation analysis

Usually, using external stimuli is a common method to explore the firing patterns of neurons. Meanwhile, we also aim to study the influence of electromagnetic induction on the firing rhythm of neuron model, so we choose $I_{ext}$ and $k_0$ as bifurcation parameters. Then Codimension-one and ISIs bifurcation are used to explore the bifurcation structure of the system (1).

3.1.   Codimension-one bifurcation analysis

The dynamical properties of external forcing current are first investigated. Equilibrium is an important part in understanding the behavior of a dynamical system. To find the equilibrium, we can rewrite the system (1) as following:

where

where $m_\infty$, $s_\infty$, $n_\infty$, $\rho\left(\varphi\right)$ are the same as defined in Section 2.

Let $F_1 = F_2 = F_3 = F_4 = 0$, then we get the expressions of $n$, $c$, $\varphi$ with respect to $V:\ n\left(V\right) = n_\infty\left(V\right)$, $c\left(V\right) = {\theta g_{Ca}{m_\infty\left(V\right)}^2\left(V_{Ca}-V\right)}/{k_{PMCA}}$, $\varphi\left(V\right) = {k_1V}/{k_2}$. Next, we put these into $F_1 = 0$, so we can easily obtain the expression of equilibrium with respect to $I_{ext}$:

We can calculate the codimension-one bifurcation points by using MATCONT, see Figure.1(a). It can be seen there are three bifurcation points on the curve, namely H, ${\rm {LP}}_1$ and ${\rm {LP}}_2$. H represents Hopf bifurcation, where $I_{ext} = -0.196411$, ${\rm {LP}}_1$ and ${\rm {LP}}_2$ are fold bifurcation points at $I_{ext} = 0.831046$, $I_{ext} = 0.703546$, respectively. The blue curve indicates stable fixed points, and the black curve is unstable fixed points. For $I_{ext}<-0.196411$, the system (2) has only one stable equilibrium and there exist quiescent state; Then the system (2) changes its stability from stable to unstable through Hopf bifurcation at $I_{ext} = -0.196411$, it means the system (2) enter firing state; For $0.703546<I_{ext}<0.831046$, the system (2) has three equilibriums through fold bifurcation ${\rm {LP}}_2$, then another fold bifurcation occurs at ${\rm {LP}}_1$ as we increase the value of parameter $I_{ext}$; Finally, for $I_{ext}>0.831046$, there is only one stable equilibrium.

Now, we only consider the effect of intensity of induced current and set $I_{ext}$ to zero to obtain the S-type bifurcation curve with respect to parameter $k_0$ in the same way, as shown in Figure.1(b). there are two fold bifurcation points on the curve occurring on $k_0 = -0.007763$ and $k_0 = 2.54\times{10}^{-22}$(not drawn). Hopf bifurcation point appears at $k_0 = 0.012850$. The bifurcation structure is similar to the previous ones, so it will not be explained in detail.

The transformation between the quiescent state and the firing state is related to these bifurcation points. Therefore, we need to verify the dynamical properties of bifurcation points. The first Lyapunov coefficient determines whether Hopf bifurcation is supercritical or subcritical. It is clear that Hopf bifurcation is supercritical (or subcritical), if the first Lyapunov coefficient is negative (or positive). We take parameter $I_{ext}$ as an example to verify the dynamical properties of these bifurcation points.

Typically, we need to calculate the Jacobian matrix of the system (2), and the Jacobian matrix can be express as:

For the bifurcation points, we use the matrix (4) and the system (2) to calculate the equilibrium, then we can get the coordinate of bifurcation points. When $I_{ext} = -0.196411$, the coordinate of equilibrium of the system (2) is $\left(V,n,v,\varphi\right) =$(-39.709558, 0.006939, 0.98244, -13.236519), hence we can obtain the corresponding eigenvalues of the Jacobian matrix $A$ are $\lambda_1 = -17.850777$, $\lambda_2 = -2.785812$, and $\lambda_{3,4} = \pm0.752697i$. Therefore, Hopf bifurcation occurs at point H. For $I_{ext} = 0.831046$, the coordinate of equilibrium of the system (2) is (-46.262568, 0.00307, 0.454643, -15.420856). And the corresponding eigenvalues of the Jacobian matrix $A$ are $\ lambda_1 = -29.52063$, $\lambda_2 = -2.96421$, $\lambda_3 = 11.504217$ and $\lambda_4 = 0$. Consequently, there is a fold bifurcation. So ${\rm {LP}}_1$ is tested. At point ${\rm {LP}}_2$ (when $I_{ext} = 0.703546$), the equilibrium of the system (2) is (-59.325527, 0.000601, 0.078164, -19.775176). the corresponding eigenvalues of the Jacobian matrix $A$ are $\lambda_1 = -33.030555$, $\lambda_2 = -2.7488253$, $\lambda_3 = 0.00146518$ and $\lambda_4 = 0$. Accordingly, a fold bifurcation occurs at ${\rm {LP}}_2$.

Besides, it is necessary to calculate the first Lyapunov coefficient [10], since it determines the stability of Hopf bifurcation. First, we can calculate the Jacobian matrix at point H:

which has a pair of conjugate eigenvalues $\lambda_{3,4} = \pm i\omega$, where $\omega = 0.752697$. Let

which satisfy $Aq = i\omega q$, $A^Tp = -i\omega p$ and $\langle p, q\rangle = 1$. Then, we need to move the equilibrium of the system (2) to the origin of the coordinate for computing the first Lyapunov coefficient. Then, we will make the transformation:

where $V_0$, $n_0$, $c_0$ and $\varphi_0$ are coordinate of H point. By this transformation, the system (2) is converted into this follow:

Now, consider the system:

where $A = \left.A\right|_H$, $F(x) = \frac{1}{2} {\boldsymbol{B}}(x, x)+\frac{1}{6} {\boldsymbol{C}}(x, x, x)+{\boldsymbol{O}}\left(\left\|x^{4}\right\|\right)$, at the neighborhood of $x$, ${\boldsymbol{B}}(x, x)$ and ${\boldsymbol{C}}(x, x, x)$ are symmetric multilinear vector functions which take on the planar vectors $x = \left(x_1,\ \ x_2,\ \ x_3,\ \ x_4\right)^T$, $y = \left(y_1,\ \ y_2,\ \ y_3,\ \ y_4\right)^T$ and $z = \left(z_1,\ \ z_2,\ \ z_3,\ \ z_4\right)^T$. In coordinates, we have

where $\xi = \left(\xi_1,\ \ \xi_2,\ \ \xi_3,\ \ \xi_4\right)$.

It is not complicated to compute

It is essential to calculate the norms in the center manifold when computing the first Lyapunov coefficient for high-dimensional. There is an invariant expression of the first Lyapunov coefficient:

Hence H is supercritical Hopf bifurcation, and it branches out the stable limit cycle. Similarly, ${\rm {H}}_1$ is also a supercritical Hopf bifurcation.

3.2.   ISIs bifurcation

The more detailed bifurcation caused by direct the current stimulation and $k_0$ can be described by ISIs series, see Figure.2. As $I_{ext}$ increases in Figure.2(a), the system (2) enters the firing state around $I_{ext} = -0.1964$. Finally, the system (2) terminates the firing state at $I_{ext}\approx0.7035$. It indicates that Hopf bifurcation and fold bifurcation are the keys to the transition between the quiescent state and the firing state.

Moreover, the system presents abundant dynamical properties for $-0.1964<I_{ext}<0.7035$. Note that in the Figure.2(a), there are complex dynamic regions around $I_{ext} = 0.0866$, $I_{ext} = 0.2636$ and $I_{ext} = 0.4956$. For $I_{ext}<0.0866$, there is the period-adding bifurcation. But the period-adding bifurcation terminates, and the maximum value of ISIs changes from increase to decrease around $I_{ext} = 0.0866$. Chaotic bursting is transformed into chaotic spike by inverse period-double bifurcation, and then turns into spike. Then an inverse period-double bifurcation occurs around $I_{ext} = 0.2636$ and leads to chaos. Subsequently, the occurrence of a shrinkage [6] causes the change of bifurcation structure and chaos disappears, which might be caused by a crisis [24]. Similarly, another period-doubling leading to chaos around $I_{ext} = 0.4956$. Here, spike is transformed into chaotic spike caused by period-double bifurcation, and finally become bursting.

The Lyapunov exponents of the dynamical system is one of the several indicators used to judge whether the dynamical system is chaotic. If the largest Lyapunov exponent is positive and the second Lyapunov exponents is equal to zero, then the system is chaotic. Accordingly, we calculate the first and second Lyapunov exponents to illustrate the chaotic behavior, as shown in Figure.3. Obviously, it is consistent with the behaviors we have described above. Note that, there is a situation that is different from the others. That is, when $I_{ext}>0.524$, the first Lyapunov exponent is positive, but the second Lyapunov exponent is negative. It can't say this situation is chaotic, since the second Lyapunov exponent is not equal to zero.

The bifurcation diagram of parameter $k_0$ in Figure.2(b) is similar to the part of bifurcation diagram of parameter $I_{ext}$. They can generate the same firing modes. It indicates that the feedback current of electromagnetic induction has the same effect on the neuron model as the external forcing currents, and the induced current acts as the external forcing current.

4.   Analysis of firing patterns

In this section, we discuss the firing patterns caused by external forcing current and $k_0$. Note that, $\rho\left(\varphi\right)$ the coupling relation between membrane potential of the neuron and magnetic flux, then we discuss the effects of strong coupling strength and weak coupling strength on the firing activities of the neuron model.

4.1.   The effect of weak coupling strength

When $\beta = 0.0001$, the feedback of magnetic flux is relatively weak on membrane potential. In this case, there are four bursting patterns, as shown in Figure.4. As can be seen from the Figure.4, the firing patterns caused by the induced current are the same as the firing patterns generated by the external forcing current. In both cases, the neuron model shows regular bursting and spike. We analyze the bursting patterns by using fast-slow dynamics analysis. Since the slow variable $c$ is distinctly slower than the other three variables in the system (2), then we can obtain the fast subsystem:

and the slow subsystem:

Thus, we treat slow variable c as a bifurcation parameter of the fast subsystem (6), then, we get four bursting types according to the bifurcation of the transition between the quiescent and firing state through numerical calculations. We present the numerical results of bifurcation analysis of the fast subsystem (6) by using MATCONT software.

4.1.1.   "Fold/fold" bursting via the "point-point" hysteresis loop

For $I_{ext} = -0.1$, we draw the bifurcation diagram of the fast subsystem with respect to the slow variable, as shown in Figure.5(a). The bifurcation curve of the fast subsystem presents Z-shape. The stable equilibria lose stability via ${\rm {H}}_1$, , and branch a stable limit cycle. Then the unstable equilibria transform into stable via ${\rm {H}}_2$, and the stable limit cycle vanishes. Note that, with the decrease of the slow variable, the lower quiescence state disappears via fold bifurcation ${\rm {LP}}_2$ on the invariant circle (SNIC), and the trajectory doesn't pass through the stable limit cycle corresponding to the firing state. After that, the trajectory attenuation oscillates to another fold bifurcation ${\rm {LP}}_1$, and then jumps into the lower quiescent state. Thus, the bursting process is completed and starts another. Since the trajectory doesn't cross the limit cycle attractor, then we only consider the bifurcation generated by the hysteresis loop. According to the classification of bursting types [8], this bursting type belongs to the "fold/fold" hysteresis loop bursting of point-point type.

4.1.2.   "Fold/Hopf" bursting via the "fold/fold" hysteresis loop

When $I_{ext} = -0.03$, the second bursting pattern is shown in Figure.5(b). The bifurcation curve of the fast subsystem with respect to slow variable $c$ is also a Z-shape curve, and the upper and lower branch of Z-shape curve intersects with the middle branch forming the fold bifurcation point ${\rm {LP}}_i\ \left(i = 1,2\right)$. Equally, there are two Hopf bifurcation points (The left-hand is not drawn) and the limit cycle is located between Hopf bifurcation points. Now we analyze the generating process of bursting through Figure.5(b). The quiescent state vanishes via fold bifurcation ${\rm {LP}}_2$, and switches to the firing state corresponding to the limit cycle attractor. Then $c$ begins to decrease and the trajectory crosses the Hopf bifurcation point ${\rm {H}}_1$ indicating the system ends the firing state. Finally, the trajectory jumps into the lower quiescent state via ${\rm {LP}}_2$. Hence the bifurcation of transformation between quiescent and firing state is "fold/Hopf" type bursting, and the bifurcation of the hysteresis loop is "fold/fold" type. Therefore, the model presents "fold/Hopf" bursting via the "fold/fold" hysteresis loop.

4.1.3.   "Fold/homoclinic" bursting via the "fold/homoclinic" hysteresis loop

Now we analyze the third type of bursting of the system (2). There are regular and chaotic bursting of this type for different $I_{ext}$.

When $I_{ext} = 0.21$, the model exhibits the "fold/homoclinic" regular bursting. As shown in Figure.5(c), the limit cycle of the upper branch of the Z-shape bifurcation curve intersects with the middle branch of the Z-shape curve forming a saddle homoclinic orbit. According to the trajectory superimposed on Figure.5(c), the model switches to the firing state through fold bifurcation ${\rm {LP}}_2$, then the repetitive spiking oscillates near to the saddle homoclinic orbit, so the firing state disappears via saddle homoclinic bifurcation. We can see the same bifurcation form the hysteresis loop. Therefore, the model exhibits "fold/homoclinic" bursting via the "fold/homoclinic" hysteresis loop.

Moreover, we also obtain the chaotic bursting of the model for $I_{ext} = 0.2406$, where the largest Lyapunov exponent is positive, and the chaotic behavior is shown by fast-slow dynamics analysis in Figure.5(d). We find that the model enters the repetitive spiking state, when the trajectory is closed to saddle homoclinic bifurcation. After a while, the model terminates the spiking state via the saddle homoclinic orbit. So, the bifurcation of transformation between the firing and quiescent and hysteresis loop is the same as $I_{ext} = 0.21$. Thus, it is "fold/homoclinic" chaotic bursting via the "fold/homoclinic" hysteresis loop.

4.1.4.   "SubHopf/homoclinic" bursting via the "fold/homoclinic" hysteresis loop

Interestingly, by changing parameter, the model can also exhibit the intrinsic pseudo-plateau bursting, which is very common in endocrine model, as shown in Figure.6.

In Figure.6(b), the quiescent state vanishes via fold bifurcation ${\rm {LP}}_2$. Then, the trajectory starts to oscillate around the stable equilibrium on the upper branch, soon afterward the trajectory crosses the subcritical Hopf bifurcation point entering the spiking state. At last, the spiking state disappears through a homoclinic orbit and shifts to the lower branch corresponding to the quiescence state. Thus, this electrical activity of the model is called "subHopf/homoclinic" bursting via the "fold/homoclinic" hysteresis loop.

4.2.   The effect of strong coupling strength

Now we consider increase the feedback of the magnetic flux on the membrane potential, and set $\beta = 0.001$ and $k_2 = 2.5$. Firstly, electrical activities of the neuron model show the firing rhythm only caused by electromagnetic radiation is different from above, and numerical results are shown in Figure.7. When the induced current is small, the neuron presents spike state. With the increase of parameter $k_0$, the firing pattern of the neuron changes into a period-2 pattern, while the neuron model enters into mixed bursting when $k_0$ increase to 0.007408. For higher $k_0$, the neuron model presents subthreshold oscillations and quiescent state.

When the combined effects of induced current and the external forcing currents are considered, the electrical activities of neuron model become more abundant. As shown in Figure.10 and Figure.12, ISIs bifurcation diagram is more complex. By changing the parameters $I_{ext}$, the neuron model can select different firing modes by applying external forcing current under induced current of different intensity in Figure.9 and Figure.11. The firing activities show that system fires with subthreshold oscillations, regular bursting, and mixed bursting for $k_0 = 0.01$. While increasing the intensity of induced current, the neuron model presents subthreshold oscillations and MMOs for $k_0 = 0.02$. Note that the regular mixed-mode oscillations alternate with the irregular mixed-mode oscillations, as shown in Figure.11 (b) and (c). Therefore, with the increase of induced current, the response of the neuron model to the stimulus current changes from bursting pattern to MMOs. It indicates that the firing behaviors of the neuron model can be regulated by both induced current and direct current stimulation.

5.   Codimension-two bifurcation analysis

In this section, we use the bifurcation theory of dynamical systems and numerical simulation to investigate the codimension-two bifurcation in the improved neuron model in $\left(I_{ext},k_0\right)$ phase plane and other parameters are given in Table 1.

5.1.   Bifurcation in $\left({\boldsymbol{I}}_{{\boldsymbol{e x t}}}, {\boldsymbol{k}}_{0}\right)$-plane

The existence of stable limit cycles indicates that the neuron model will present firing activities, such as spike or bursting. The region of oscillation activities can be detected by two-parameter bifurcation analysis, as shown in Figure.13. The region of oscillation activities roughly locates in the region surrounded by Hopf bifurcation curve $h_1$ in Figure.13(b). For fixed $k_0$, the bifurcation of the transformation between the quiescent state and the firing state is related to the boundaries of the region surrounded by $h_1$. For example, when $k_0 = 0.03$, the boundary points are Hopf bifurcation points. It means that the bifurcation of the transformation between the quiescent state and the firing state is Hopf bifurcation. While $k_0 = 0.01$, the bifurcation is fold bifurcation and Hopf bifurcation.

We use MATCONT software to calculate codimension-two bifurcation points. The meaning of each label in the figure is interpreted as follows: ${\rm {CP}}_i\left(i = 1,\ldots, 6\right)$ are the cusp bifurcation; ${\rm{ BT}}_i(i = 1,2,3)$ represent Bogdanov-Takens bifurcation; ${\rm {ZH}}_i\left(i = 1,\ 2\right)$ indicate fold-Hopf bifurcation; ${\rm {GH}}_i\left(i = 1,\ldots,10\right)$ represent Bautin bifurcation; HH is Hopf-Hopf bifurcation; ${\rm {NS}}_i\ \left(i = 1,\ 2,\ 3\right)$ are neutral saddle. Readers can refer to [10] for understanding the meaning of this bifurcation points. The more detailed data is shown in Table 2.

From Figure.13 (a) and (b), we can see that the fold bifurcation curve $f_1$ and Hopf bifurcation curve $h_1$ seem to be coincident into one curve, and two curves are tangent at ${\rm{ BT}}_2$. It indicates that the equilibrium and eigenvalue change greatly, even if the parameter value changes slightly. The same situation occurs between the curves $f_2$ and $h_2$, and they are tangent at ${\rm{ BT}}_3$. Moreover, most of the codimension-two bifurcation points occur on $h_1$ and $f_2$.

For cusp points ${\rm {CP}}_i$, the system (2) has a zero eigenvalue. Especially, near each ${\rm {CP}}_i$, the system (2) is locally topologically equivalent to the following normal form:

where $\xi_{\pm} \in {\boldsymbol{R}}^{1}, {\boldsymbol{R}}^{2}, \beta_{1,2} \in \bf{R}$

There are three Bogdanov-Takens bifurcations labeled ${\rm{BT}}_i\ \left(i = 1,\ 2,\ 3\right)$ with two eigenvalues $\lambda_{1,2} = 0$. Bogdanov-Takens bifurcation is the tangency point of the Hopf bifurcation curve and fold bifurcation curve. Near the point ${\rm{BT}}_i$, the system can be simplified as the following topological normal form:

where $\beta_1,\ \beta_2\in{\boldsymbol{R}}$, $s = sign(ab) = 1$.

There are ten generalized Hopf bifurcation points. We know Bautin bifurcation is the degenerated Hopf bifurcation whose first Lyapunov coefficient is equal to zero. Near the Bautin bifurcation points ${\rm {GH}}_i\left(i = 1,\ldots4,\ 6,\ldots,10\right)$, the differential system is equivalent to the following topological normal form:

where $\beta_1,\beta_2\in{\boldsymbol{R}}$, when real eigenvalue is positive, ${\dot{\xi}}_+ = \xi_+$, otherwise ${\dot{\xi}}_- = \xi_-$,

Near the bifurcation point ${\rm {GH}}_5$, the differential system is equivalent to the following topological normal form:

where $\beta_1,\ \beta_2\in{\boldsymbol{R}}$, $s = sign\ \left(l_2\right) = -1$.

Further, there are two fold-Hopf bifurcation points, their eigenvalues consist of one zero eigenvalue, a pair of purely imaginary eigenvalues, and one nonzero real eigenvalue, but they don't have the fixed normal form. A Hopf-Hopf bifurcation labeled HH occurs, and it has two pairs of purely imaginary eigenvalues, but it has no fixed normal form like the fold-Hopf bifurcation. The meaning of ${\rm {NS}}_i \left(i = 1,\ 2,\ 3\right)$ are zero-neutral saddle with one zero eigenvalue and two real eigenvalues that satisfy their sum is equal to zero.

5.2.   Bogdanov-Takens bifurcation analysis

In this section, we investigate the Bogdanov-Takens bifurcation of the system (2) through the method proposed by [3]. We treat $\left(I_{ext},k_0\right)$ as a pair of bifurcation parameters, other parameters are given by Table 1. We select ${BT}_2$ as the point of analysis, where $\left(I_{ext},k_0\right) = \left(0.649385813300958,\ \ 0.00912704116143242\right) = \mu_0$, and its coordinate is $\left(V,n,c,\varphi\right) = \left(-60.0447105,\ 0.0005496996,\ 0.07053964,\ -20.014904\ \right) = x_0$. Now we can rewrite the system (2) as the following:

where $x = \left(V,n,c,\varphi\right)$, $\mu = \left(I_{ext},k_0\right)$, and

where $m_\infty$, $s_\infty$, $n_\infty$, $\rho\left(\varphi\right)$ are the same as before.

Afterward, we can obtain the Taylor series of $F\left(x,\mu\right)$ around $\left(x_0,\mu_0\right)$ as the following:

Note

Then we can get the eigenvalues of matrix $B$, specifically, 0, 0, -33.09149, -2.765521. Let $P = \left(p_1,p_2,P_0\right)$ be an invertible matrix, which satisfies $P^{-1}AP = \left(\begin{matrix}J_0&0\\0&J_1\\\end{matrix}\right)$, where $J_0 = \left(\begin{matrix}0&1\\0&0\\\end{matrix}\right)$, $J_1 = \left(\begin{matrix}-33.0914934&0\\0&-2.765521027\\\end{matrix}\right)$, $p_1$, $p_2$ are the generalized eigenvectors of matrix $B$ corresponding to the double-zero eigenvalue and $P_0$ consists of the generalized eigenvectors of the matrix $J_1$.

So, we can get

We define $P^{-1} = \left(q_1,q_2,Q_0\right)^T$, then

According to expression (28) and (29) in [3], we can obtain

Let $\lambda_1 = I_{ext}-0.64939$ and $\lambda_2 = k_0-0.009127$, then we treat $\lambda_1$, $\lambda_2$ as a pair of bifurcation parameters. Thus, we have

By using Theorem 1 in [3], the system (12) at $x = x_0$, $\mu = \mu_0$ is locally topologically equivalent to

Moreover, we can transform the variables by

Then system (12) is transformed into the following:

where

By the theory in [10], we can compute the following terms:

According to the theory in [10] and our results of calculation, we have the following:

Theorem 5.1. Let $\left(I_{ext},k_0\right)$ be the Bogdanov-Takens bifurcation point of systems (12). Denote ${\bar{\lambda}}_1 = I_{ext}-0.64939$, ${\bar{\lambda}}_2 = k_0-0.009127$. Then if $0<{\|\bar{\lambda}}_1,{\bar{\lambda}_2\|}^2\ll1$, the dynamics on the center manifold of the system (12) near $x = x_0$, $\mu\approx\mu_0$ is locally topologically equivalent to the following system:

In addition, system (16) has three bifurcation curves in the small neighborhood near the origin:

(i)there is a saddle-node bifurcation curve:

(ii) there exists a Hopf bifurcation curve:

(iii) there exists a saddle homoclinic bifurcation curve:

6.   Conclusion

The detailed bifurcation information on the improved endocrine model is presented under electromagnetic induction and external current in this paper. For external current and parameter $k_0$, we use codimension-one and ISIs bifurcation to study the bifurcation structure. Numerical simulations show that the system realizes the transition between quiescence and firing state via Hopf bifurcation and fold bifurcation. As the $I_{ext}$ increases, the system experiences three complex dynamical behaviors: period-adding bifurcation, inverse period-doubling bifurcation, and period-doubling bifurcation leading to chaos. It reveals that external forcing current may have significant influences on the improved endocrine model under the magnetic field, and the induced current plays the role of external forcing current. The numerical results show that the same discharge modes can be generated by induced current and external forced current when the coupling strength between magnetic flux and membrane potential is weak. Besides, we also analyze the different bursting patterns by fast-slow dynamics analysis. While the coupling strength between magnetic flux and membrane potential is strong, the neuron model exhibits the different types of electrical activities under different intensity of induced current. This indicates that the electrical activities of the model neuron are very sensitive to the response of induced current and stimulus current. Beyond that, we explore the region of oscillation activities by using two-parameter bifurcations and calculate the local topologically equivalent normal forms. If the fold bifurcation and Hopf bifurcation curve are too close, it is not easy to draw the homoclinic bifurcation curve. Nevertheless, we also discuss the dynamical properties near the Bogdanov-Takens bifurcation point in detail, and theoretically confirm the three bifurcation curves by using bifurcation theory and central manifold theorem. Our theoretical and numerical results may contribute to better understand the dynamical behavior and the mechanism of the endocrine neuron under the magnetic field and external current terms.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant Nos. 11872183 and 11572127.