Influence of Electric Current and Magnetic Flow on Firing Patterns of Pre-Bötzinger Complex Model

The dynamics of neuronal firing activity is vital for understanding the pathological respiratory rhythm. Studies on electrophysiology show that the magnetic flow is an essential factor that modulates the firing activities of neurons. By adding the magnetic flow to Butera's neuron model, we investigate how the electric current and magnetic flow influence neuronal activities under certain parametric restrictions. Using fast-slow decomposition and bifurcation analysis, we show that the variation of external electric current and magnetic flow leads to the change of the bistable structure of the system and hence results in the switch of neuronal firing pattern from one type to another.


Introduction
Breathing is an important physiological activity that is necessary for all mammals, including human beings, to sustain their lives. Experiments have shown that respiratory rhythms in the neonatal nervous system of mammals may be related to pacemaker neurons in the pre-Bötzinger complex (pre-BötC) [1,2]. Solomon et al. found that the chemical stimulation of the pre-BötC in vivo manifests respiratory modulation consistent with a respiratory rhythm generator [3]. The phrenic motor activity evoked by chemical stimulation of the pre-BötC is affected and integrated through the regulation of the respiratory network driven by input from central and peripheral chemoreceptors [4]. Besides, Koizumi and Smith used real-time calcium activity imaging combined with whole-cell patchclamp recording to analyze contributions of subthreshold conductances in the excitatory rhythm-generating network [5]. Negro et al. demonstrated that dendritic Ca 2+ activates an inward current to electronically depolarize the soma, rather than propagating as a regenerative Ca 2+ wave [6]. Research evidences underscored that respiratory rhythmogenesis may depend on dendritic burstgenerating conductance activated in the context of network activity.
The mathematical framework of neuron electrophysiological models has been derived from the Hodgkin-Huxley (H-H) model established by Hodgkin and Huxley [7]. Based on the H-H model, Butera et al. created two mathematical models of inspiratory pacemaker neurons [8,9] to simulate the respiratory rhythm generation of single as well as coupled pre-BötC neuron. Subsequent studies have shown that neurons within the pre-BötC have a persistent sodium (NaP) current and a calcium-activated nonspecific cationic (CAN) current. CAN current can be activated via secondmessenger-mediated synaptic pathways [10,11]. Driven by the experimental results, Toporikova and Butera proposed a pre-BötC (TB) model that encompasses both I NaP -dependent somatic bursting and I CAN -dependent dendritic bursting [12]. The model explains a number of conflicting experimental results, and it is able to generate a robust bursting rhythm, over a large range of parameters, with a frequency adjusted by neuromodulators. Later, Park and Rubin simplified the TB model and proposed a singlecompartment model of a pre-BötC inspiratory neuron which is able to generate all major activity patterns seen in the twocompartment model [13]. Also, the neuronal models exhibit abundant dynamic characters such as the bifurcation phenomenon. Gu et al. demonstrated bifurcation and complex dynamic behaviors in biological experiments [14,15].
An appropriate external stimulus can change the firing patterns of neurons. Much research in recent years has focused on the effects of electromagnetic radiation on neuronal behaviors [16][17][18][19][20]. Electromagnetic radiation can affect the dynamic characteristics of neurons, and electrical or electromagnetic stimulation can also be used to treat neurological diseases [21][22][23]. Recently, Song et al. designed a nonlinear circuit including an inductor, resistor, capacitor, and other electric devices. They found that the energy storage is dependent on the external forcing and the energy release is associated with the electric mode [24]. Lv and Ma also found that electromagnetic radiation can not only excite quiescent neurons but also suppress the electrical activities in the improved three-variable Hindmarsh-Rose model [25]. Duan et al. added magnetic flow as a new variable to the Butera model to explore the effects of electromagnetic induction on neuronal activities [26]. Parastesh et al. proposed a new memory function based on discontinuous flux coupling and studied various dynamic characteristics of discontinuous flux coupling neuron models [27]. Considering the mutual influence of electric and magnetic fields, Rostami and Jafari described a new Hindmarsh-Rose (HR) neuron model with more bifurcation parameters to study the formation of defects in excitable tissues and the resulting emission waves [28]. Although there have been some studies on the influence of electromagnetic fields on neuronal activity patterns, little attention has been paid to the influence of electromagnetic currents on pre-BötC.
This paper considers the effect of electromagnetic induction on the Butera neuron model with external electric current, magnetic flux, and CAN current. The membrane potential can be adjusted by a memristor that connects the membrane potential to the magnetic flux [29][30][31][32]. This paper is organised as follows. In Section 2, we describe the Butera model with a memristor. Based on this model, we explore the firing patterns of the system in Section 3. The effects of external forcing current and the magnetic flux on the bursting rhythms of the system are studied by nondimensionalization analysis, fast-slow decomposition methods, and twoparameter bifurcation analysis. The dynamical mechanisms of generation and transition of firing patterns are given. The conclusion is given in the last section.

Model Description
Electric and magnetic flow is introduced in Butera's singlecompartment model of pre-BötC. The model is described as follows: where V is the membrane potential and n and h are gating variables for the voltage-gated potassium and sodium channels, respectively. The functions I Na , I K , I L , I NaP , I CAN , and I tonic−e represent fast sodium current, delayed rectifier potassium current, leakage current, persistent sodium current, calcium-activated nonspecific cationic current, and external tonic drive current, respectively. I extz is a direct current of external stimuli. C is the whole cell capacitance. Particularly, The calcium dynamics is given as [13] d Ca where l represents the fraction of IP 3 channels in the membrane of endoplasmic reticulum (ER) that have not been inactivated, which depends on the intracellular calcium concentration ([Ca]). Equation (3a) specifies that [Ca] is determined by the flux into the cytosol from the ER (J ER IN ) and the flux out of the cytosol into the ER (J ER OUT ). These fluxes are regulated by the intracellular concentration of IP 3 , [IP 3 ], and IP 3 channel gating variable, l. The variable φ refers to the magnetic flux across the cell membrane. ρðφÞ is the coupling strength between membrane potential of neuron and magnetic flux, which is a magnetic flux-controlled memristor, and it is equivalent to memory conductance: ρðφÞ = α + 3βφ 2 , where α and β are fixed parameters [29,30]. k 1 and k 2 show the relationship between the membrane potential and the magnetic flux. And we suppose that I app = I extz − k 1 VρðφÞ in Equation (1a). The term k 1 VρðφÞ introduces the inhibitory modulation of membrane potential as induction current results from variation of magnetic flux and field [25,33], and it can be described as follows [33]: Dynamical analysis indicates that subsystem (3a) and (3b) is in an active state only when IP 3 lies between 0.95 μM and 1.4 μM [13]. Therefore, in this study, we set ½IP 3 = 1:2 μM, which implies that ½Ca is not invariant. The function expressions and parameter values of other variables are presented in Appendix.

Main Results
We used the method of fast-slow decomposition to study the firing patterns of system (1a)-(3b). In order to clearly identify the timescales of different variables, we nondimensionalize the full system (1a)-(3b) as that done in previous work [34].

Nondimensionalization and Simplification of Timescales.
The variables are rescaled so as to further reveal their timescales. To this end, we define new dimensionless variables (v, c, φ, τ) and voltage, calcium, magnetic flux, and time scales Q v , Q c , Q φ , and Q t , respectively, such that Note that n, h, and l are already dimensionless in Equations (1b), (1c), and (3b).
Since the membrane potential V typically lies between −60 mV and 0 mV, we define T x = max ð1/τ x ðVÞÞ over the range V ∈ ½−60, 0 and then define t x ðVÞ, a rescaled version of τ x ðVÞ, by t x ðVÞ = T x τ x ðVÞ for x ∈ fn, hg. As ½Ca typically lies between 0 μM and 1 μM, we define Gð½CaÞ = ½IP 3 ½Ca/½ð½IP 3 + K I Þð½Ca + K a Þ, g SERCA ð½CaÞ = V SERCA ð½Ca/K 2 SERCA + ½Ca 2 Þ over the range ½Ca ∈ ½0, 1. Furthermore, we set G c = max ðG 3 ð½CaÞÞ and G S = max ðg SERCA ð½ CaÞÞ and then define P max = max fL IP 3 , P IP 3 G c , G S g. We also define g max = max g Na , g K , g L , g NaP , g tonic−e , g CAN , I extz , k 1 f g : ð6Þ According to system (1a)-(3b), we get the following dimensionless system: Combining the ranges of V, φ, and ½Ca, suitable choices for the voltage, magnetic flux, and calcium scales are Q v = Q φ = 100 mV and Q c = 1 μM. We also see that values of m ∞ ðVÞ, mp ∞ ðVÞ, n ∞ ðVÞ, h ∞ ðVÞ, f ð½CaÞ, GðCaÞ, g SERCA ð½ CaÞ, n, h, and l all lie in the range ½0, 1. Combining the parameter values in Table 1 and the values of variables in the paper, we have g max = max fI extz g = 50 μA. Specifically, Figures 1(a) and 1(b) show the plots of 1/τ n ðVÞ and 1/τ h ðV Þ over the range V ∈ ½−60, 0, which indicates that T n ≈ 1:3 ms −1 and T h ≈ 0:0025 ms −1 . Similarly, we obtain G c ≈ 0:06 and G S ≈ 1000 pL · ms −1 from Figures 1(c) and 1(d), respectively. So, we have P max ≈ 1860 pL · ms −1 . Using these values, we see that all terms in the right-hand sides of Equations (7a)-(7f) are bounded (in absolute value) by one.
The coefficients of the derivatives on the left sides of Equations (7a)-(7f) now reveal the relative rates of evolution of the variables. We find that C/g max = 0:42 ms~Oð1Þ ms, 1 /T n = 0:77 ms~Oð1Þ ms, 1/T h = 400 ms~Oð100Þ ms, σ/ð P max · K Ca Þ ≈ 3:98ms~Oð10Þms, and 1/Q c · A = 200 ms~Oð 100Þ ms. We use the notation O to denote an order of magnitude estimate: x~ð10 n Þ, where n is the nearest integer to log ðxÞ. We choose the fast timescale as our reference time, i.e., pick Q t = 1 ms, and set As a result, dimensionless system (7) becomes system (9), Neural Plasticity namely, with relative rates of all variables: From this, we conclude that v, n, and φ evolve on a fast timescale, c evolves on a slow time scale, and h and l evolve on a super slow time scale. Thus, the above model can be regarded as a dynamic model containing three timescales: fast, slow, and super slow variables. Equations (1a), (1b), and (1d) form the fast subsystems, Equation (3a) is the slow subsystem, and Equations (1c) and (3b) form the super slow subsystem. In the current model, intracellular calcium dynamics is confined to subsystem (3a) and (3b) and evolves independently from subsystem (1a)-(1d). The dynamics of [Ca] is affected by intracellular Ca 2+ concentration [Ca] and IP 3 channel gating variable l (Equations (3a) and (3b)). For simplicity, we chose g CAN Tot = g CAN f ð½CaÞ as the bifurcation parameter, where f ð½CaÞ is a monotonically increasing concave function of [Ca].
According to the nondimensionalization of the model, the variable h is super slow. So, h can be considered a constant when we do the fast-slow decomposition. That is, we take h as the average value of the variable in the following analysis.

The Firing Patterns and Bifurcation Analysis without
External Stimulation. Firstly, we consider the firing activity of the full system with zero electric current and magnetic flow. Time courses of the membrane potential V (black solid) and the intracellular calcium concentration [Ca] (red dashed) are shown in Figure 2(a).
One-parameter and two-parameter bifurcation diagrams of the fast subsystem are shown in Figures 2(b) and 2(c), respectively. According to the nondimensionalization, the variable h is a super slow variable. So, h is con-sidered to be a constant which we take its average value as h = 0:2834. g CAN Tot is a slow variable and regarded as a bifurcation parameter of the fast system. In this figure, the equilibrium points form an S-shaped curve. The lower (black solid line) and middle (black dash-dot line) branches of the curve are composed of the stable nodes and unstable saddles, respectively. The upper branch of the curve is composed of the stable and unstable focus separated by the subcritical Hopf bifurcation point (subH). With the decrease of g CAN Tot , the stable focus becomes unstable, and an unstable limit cycle (red dashed line) occurs at the subcritical Hopf bifurcation (subH). The unstable limit cycle generated by the subcritical Hopf bifurcation transits into the stable limit cycle (red solid line) through the fold bifurcation of limit cycle (LPC). The points F 1 and F 2 refer to the fold bifurcation. The stable limit cycle disappears due to the homoclinic bifurcation (HC). The trajectory of the full system (green curve) is also appended on the bifurcation diagram. According to Izhikevich's classification scheme of bursting [35], the firing pattern of system (1a)-(3b) can be identified as "subHopf/homoclinic" bursting via "fold/homoclinic" hysteresis loop. Figure 2(c) shows the two-parameter bifurcation analysis of the fast subsystem in (g CAN Tot ,h)-plane, where fc, hc, lc, and homo represent the fold bifurcation curve (red), Hopf bifurcation curve (blue), fold limit cycle bifurcation curve (black), and homoclinic bifurcation curve (purple), respectively. The trajectory of the full system curve (green) is also appended in the (g CAN Tot ,h)-plane. With the increase of g CAN Tot , the trajectory passes through different bifurcation curves. Then, with the decrease of g CAN Tot , the trajectory crosses same bifurcation curves again and drops to the origin state after meeting the homoclinic curve.
3.3. The Effect of Electric Current on Activity of Pre-BötC. Figure 3 shows some electrical activities of the system with different external forcing current I extz . For other parameter values, please refer to the appendix.
The enlarged area of the two periods of Figure 3(a) is shown in Figure 4(a). We divide the change of membrane potential with time in a cycle into phases ①-④. Phase ① (from ★ to ▲) is the low potential resting stage, and its duration is about 1200 ms. Phase ② (from ▲ to ■) is tonic spikings with gradually decreasing amplitude. The time interval between adjacent spiking is very short, close to zero. Phase ③ (from ■ to ◆) is spikings with small amplitude, and the interval between spiking is very short. Phase ④ (from ◆ to ★) is spikings with gradually increasing amplitude, and the interval between adjacent spiking is also very short.
In order to further reveal the influence of external forced current on the electrical activity of neurons, the ISI (interspike interval) bifurcation diagram is drawn in Figure 4(b). We take the time interval between the maximum values of adjacent spiking as an ISI bifurcation diagram and use the logarithmic function to process the data. While the duration of phase ① increases slowly with the  Neural Plasticity increases of current I extz at first, it begins to decrease rapidly when the current value reaches 21.5 μA/cm 2 . At the same time, the start time of phase ② advances, which leads to more spiking in the original low potential resting state, as shown in Figure 2 from (c) to (d). When the current increases to 33.5 μA/cm 2 , the original low potential resting state disappears and all become spiking. In the process of current I extz increase, phase ③ maintains small amplitude spikings and the duration of time gradually decreases. When the value of I extz reaches 47.7 μA/cm 2 , the original spiking of phase ③ becomes a resting state with a higher membrane potential, and then, the duration gradually increases, resulting in a gradual increase in the current value near 50 μA/cm 2 of an ISI sequence.

Fast-Slow Decomposition.
We can obtain oneparameter bifurcation diagrams of the fast subsystem (1a), (1b), and (1d) as shown in Figure 5, in which the projection of trajectory (the green curve) of the full system is also superposed. The equilibrium points form an Sshaped curve which is similar to that without the electric current and magnetic flow. The lower (solid black line) and middle (dash-dot line) branches of the curve are composed of stable nodes and unstable saddles, respectively. The upper branch of the curve is composed of stable and unstable focuses separated by the subcritical Hopf bifurcation point (subH). With the decrease of g CAN Tot , the stable focuses become unstable, and the unstable limit cycles (red dashed line) occur at the subcritical Hopf bifurcation (subH) on the upper branch. The unstable limit cycles turn to be stable limit cycles (red solid circle) via fold bifurcation of limit cycle (LPC). The points F 1 and F 2 refer to fold bifurcation, and the point HC represents homoclinic bifurcation.  For I extz = −2 μA/cm 2 and I extz = 5 μA/cm 2 , the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable g CAN Tot is shown in Figures 5(a) and 5(b), respectively. Here, we set h = 0:2788 and h = 0:2375, respectively (h takes the average value). Firstly, the trajectory transits from the lower rest state to the spiking state via fold bifurcation (F 1 ). Due to the attraction of the stable focus, the amplitude of oscillation rapidly decreases until the trajectory passes through the subcritical Hopf bifurcation (subH).
Then, as a result of the repelling effect of unstable focus, the amplitude of the trajectory increases gradually. Finally, the trajectory transits from spiking state to lower rest state, which completes one periodic oscillation. Thus, the firing activity of system (1a)-(3b) is the "subHopf/subHopf" bursting via "fold/homoclinic" hysteresis loop. In this case, the bistable structure is composed of the stable node on the lower branch and the stable focus on the upper branch of the equilibrium curve.

Neural Plasticity
When I extz increases to 25 μA/cm 2 and 30 μA/cm 2 , the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable g CAN Tot is shown in Figures 5(c) and 5(d), respectively. In these two cases, the parameter h is set to be h = 0:1288 and h = 0:0986, respectively. The trajectory jumps up to the upper branch of the S-shaped curve via fold bifurcation F 1 . With the strong attraction of stable focus, the oscillation decays rapidly around the stable focus. After passing through the subcritical Hopf bifurcation point (subH), the trajectory unfolds itself progressively due to the repelling of the unstable focus and the attracting of the stable limit cycle generated by the fold limit cycle bifurcation (LPC). Finally, the trajectory drops to the lower branch of the equilibrium curve via saddle-node bifurcation on an invariant circle. The bistable structure is composed of the stable nodes on the lower branch and stable limit cycles on the upper branch. The bursting can be classified as "subHopf/fold cycle" bursting via "fold/circle" hysteresis loop.
When the value of I extz continuously increases to 40 μA/cm 2 and 50 μA/cm 2 , the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable g CAN Tot is shown in Figures 5(e) and 5(f), respectively. h is set to be h = 0:0490 and h = 0:0307, respectively. The bistable structure in this case is confirmed by the stable limit cycles and the stable focus on the upper branch of the equilibrium curve. The quiescent state is lost via subcritical Hopf bifurcation (subH), and the periodic limit cycle attractor corresponding to repetitive spiking disappears via fold limit cycle bifurcation (LPC). Hence, the bursting is "subHopf/fold cycle" bursting.
In Figures 5(c)-5(f), periodic orbits disappear via a saddle-node bifurcation on an invariant circle (SNIC) as g CAN Tot decreases. Take I extz = 30 μA/cm 2 ( Figure 5(d)) as an example, the periodic orbit of the (½Ca,l)-space is shown in Figures 6(a) and 6(b). Let ½Ca SNIC be the value of [Ca] at SNIC; the period of the limit cycle tends to infinite when g CAN Tot ð½Ca SNIC Þ ≈ 0:02937 at which the saddle-node bifurcation on an invariant circle (SNIC) occurs. Then, we get ½Ca SNIC ≈ 0:02942 according to g CAN Tot ð½Ca SNIC Þ = g CAN f ð ½Ca SNIC Þ.The corresponding periodic orbit in (½Ca, l )-space is shown in Figure 6(b). The point P 1 (blue square) denotes the time when the trajectory in subsystem (3a) and (3b) jumps up from the left knee to the right branch of the [Ca]-nullcline. Once leaving the left knee, the trajectory jumps to the right branch quickly. Thus, we have a rapid increase of [Ca] and g CAN Tot correspondingly, and this increase pushes the trajectory into the branch of periodic orbits. Then, the value of [Ca] slowly decreases while the trajectory in (½Ca, l)-space slowly moves along the right branch of the [Ca]-nullcline and finally jumps back to the left branch of [Ca]-nullcline. As long as [Ca] is greater than ½Ca SNIC (P 2 ), the system keeps firing. Once [Ca] falls below ½Ca SNIC , the trajectory solution stops firing and jumps down to the lower branch of the S-shaped curve ( Figure 5(d)). Next, the trajectory in (½Ca, l)-space moves up along the left branch of the [Ca]-nullcline towards its original position, and this completes one cycle of the periodic orbit in the subsystem (3a) and (3b).

Bifurcation Analysis.
To better understand the mechanism of the different firing patterns evoked by the electrical current, we present two-parameter bifurcation analysis of  Figure 6: (a) The period of limit cycle of the fast subsystem with respect to parameter g CAN Tot at I extz = 30 μA/cm 2 ; (b) the periodic orbit in subsystem (3a) and (3b). The red curve is the [Ca]-nullcline, and the green curve is l-nullcline. The black closed orbit is the solution of subsystem (3a) and (3b) projected in (½Ca,l)-space. 10 Neural Plasticity For I extz = −2 μA/cm 2 and I extz = 5 μA/cm 2 , as shown in Figures 7(a) and 7(b), the trajectory of the full system passes through the fold bifurcation curve (fc), the Hopf bifurcation curve (hc), and the fold limit cycle bifurcation curve (lc) at the same time. The fold bifurcation curve and the Hopf bifurcation curve intersect near the full system trajectory.
When I extz increases to 25 μA/cm 2 and 30 μA/cm 2 , the trajectory of the full system moves down with the decreasing of h. It still passes through the three bifurcation curves of fc, hc, and lc. But the gap distance between the fold bifurcation (fc) and the Hopf bifurcation curve (hc) increases, as shown in Figures 7(c) and 7(d).
When the value of I extz continuously increases to 40 μA/cm 2 and 50 μA/cm 2 , the trajectory of the full system continues to move down to the decreasing direction of h. It passes through two bifurcation curves of hc and lc. The trajectory no longer passes through the fold bifurcation curve, which leads to the changes of the bistable structure, as shown in Figures 7(e) and 7(f).
The relative position of bifurcation curves and the trajectory of the system projected in (g CAN Tot , h)-plane changes with the increase of I extz . That means the value of I extz affects the critical bifurcations that determine the pattern of  12 Neural Plasticity bursting. With the increase of I extz , the bistable structure related to bursting is changed, and the bursting with low potential resting (the stable node on the lower branch) transits to bursting with high potential resting (the stable focus on the upper branch).

Influence of Magnetic Flow on Activity of Pre-BötC.
We discuss the effects of magnetic flow on firing patterns of the system by switching the parameter k 1 . Firing activities of the system corresponding to different values of k 1 are shown in Figure 8. Other parameter values are given in the appendix. The firing patterns of system (1a)-(3b) transit from one bursting pattern to another one with the external magnetic flow k 1 increases. The ISI bifurcation diagram of the magnetic flow feedback coefficient k 1 is shown in Figure 9. Circled numbers are the same as that in Figure 4(a). As k 1 increases, the duration of phase ① first slowly increases and then decreases. When k 1 increases to 0.88, the original low potential resting state turns to be spiking. With the increase of k 1 , phase ③ exhibits small amplitude oscillating, and its duration gradually decreases. When the value of k 1 reaches 1.33, the small amplitude oscillating in phase ③ becomes a resting state with a higher membrane potential, and at the same time, its duration gradually increases, which results in a gradual increase of ISI sequence.
Similarly, Equations (1a), (1b), and (1d) form the fast subsystem, and [Ca] is still the slow variable. We use g CAN Tot = g CAN f ð½CaÞ, instead of [Ca], as the bifurcation parameter, and set h to be constant. The bifurcation diagrams of the fast subsystem on the (g CAN Tot ,V)-plane are shown in Figure 10.

Neural Plasticity
With the parameter k 1 increasing, the bistable structure also changes. The bursting with low potential resting (the stable node on the lower branch shown in Figures 10(a) and 10(b)) transits to bursting with high potential resting (the stable focus on the upper branch shown in (Figures 10(c) and 10(d)). In (Figures 10(b)-10(d)), periodic orbits also disappear via a saddle-node on an invariant circle (SNIC) bifurcation as g CAN Tot decreases. The analysis is similar to that in Figure 6.
Two-parameter bifurcation of the fast subsystem in (g CAN Tot ,h)-plane is shown in Figure 11. Similarly, with the increase of k 1 , the distance between the curves of fold bifurcation (hc) and Hopf bifurcation (lc) increases, which results in the changes of the bistable structure. The parameter k 1 and the electric current have the same effect on the bursting transition.
With the increase of k 1 , the middle branch of the Sshaped curve moves to the right and the lower branch to the left, as shown in Figure 12. Compared with the case when there is no magnetic flow stimulus, the increase of k 1 does not change the bifurcation structure but changes the position of the bifurcation point.
Finally, we give the ISI bifurcation diagram of α and β to illustrate the influence of the magnetic flow on firing patterns, as shown in Figures 13(a) and 13(b), respectively. The bifurcation structure of ISI is similar to that of Figures 4(b) and 9, which indicates that as the value of α or β increases, the shape of the bursting will change. When the value of α or β is large enough, the bursting will transit from low potential resting to high potential resting at different threshold values. Increasing the value of α or β can accelerate the transition from lowpotential resting state to high-potential resting state.

Conclusion
Bifurcation is one of the most important tools for understanding the generation and transition of rhythms. Changes in system parameters or external stimulations lead to the appearance of different bifurcations. Based on Butera's model with a memristor, the effects of electric current and magnetic flow on firing patterns of the system are studied by means of fast-slow decomposition and bifurcation analysis. We found that both the direct current and magnetic flux affect the rhythm of the pre-BötC neuron significantly.
The structure of two-parameter bifurcation gives us much information about the pattern of firing. The direct current I extz can affect the relative position of the bifurcation curves, which will lead to changes of the bistable structure of the fast subsystem. The stability composed of the stable nodes on the lower branch and stable focus on the upper branch changes to the stability composed of stable focus and stable limit cycle on the upper branch with the increase of I extz , which indicates that the control of direct current I extz can lead to desired firing patterns of the system. The variation of magnetic flow can result in a similar effect. The

15
Neural Plasticity results in this study may help to reveal and explain the dynamic mechanism of electromagnetic pathogenesis and may provide theoretical guidance to diagnosing the diseases with pathological respiratory rhythm.