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.

1. 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 patch-clamp recording to analyze contributions of subthreshold conductances in the excitatory rhythm-generating network [5]. Negro et al. demonstrated that dendritic Ca2+ activates an inward current to electronically depolarize the soma, rather than propagating as a regenerative Ca2+ wave [6]. Research evidences underscored that respiratory rhythmogenesis may depend on dendritic burst-generating 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 second-messenger-mediated synaptic pathways [10, 11]. Driven by the experimental results, Toporikova and Butera proposed a pre-BötC (TB) model that encompasses both INaP-dependent somatic bursting and ICAN-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 single-compartment model of a pre-BötC inspiratory neuron which is able to generate all major activity patterns seen in the two-compartment 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-20]. Electromagnetic radiation can affect the dynamic characteristics of neurons, and electrical or electromagnetic stimulation can also be used to treat neurological diseases [21-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-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 two-parameter bifurcation analysis. The dynamical mechanisms of generation and transition of firing patterns are given. The conclusion is given in the last section.

2. Model Description

Electric and magnetic flow is introduced in Butera's single-compartment 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 INa, IK, IL, INaP, ICAN, and Itonic−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. Iextz is a direct current of external stimuli. C is the whole cell capacitance. Particularly, The calcium dynamics is given as [13] where l represents the fraction of IP3 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 (JER) and the flux out of the cytosol into the ER (JER). These fluxes are regulated by the intracellular concentration of IP3, [IP3], and IP3 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]. k1 and k2 show the relationship between the membrane potential and the magnetic flux. And we suppose that Iapp = Iextz − k1Vρ(φ) in Equation (1a). The term k1Vρ(φ) 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 IP3 lies between 0.95 μM and 1.4 μM [13]. Therefore, in this study, we set [IP3] = 1.2 μM, which implies that [Ca] is not invariant. The function expressions and parameter values of other variables are presented in Appendix.

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

3.1. Nondimensionalization and Simplification of Timescales

The variables are rescaled so as to further reveal their timescales. To this end, we define new dimensionless variables () and voltage, calcium, magnetic flux, and time scales , and Q, 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 = max(1/τ(V)) over the range V ∈ [−60, 0] and then define t(V), a rescaled version of τ(V), by t(V) = Tτ(V) for x ∈ {n, h}. As [Ca] typically lies between 0 μM and 1 μM, we define G([Ca]) = [IP3][Ca]/[([IP3] + K)([Ca] + K)], gSERCA([Ca]) = VSERCA([Ca]/KSERCA2 + [Ca]2) over the range [Ca] ∈ [0, 1]. Furthermore, we set G = max(G3([Ca])) and G = max(gSERCA([Ca])) and then define Pmax = max{LIP, PIPG, G}. We also define

According to system (1a)–(3b), we get the following dimensionless system: with dimensionless parameters (x ∈ {Na, K, L, NaP, CAN, tonic − e}), (x ∈ {Na, K, L, syn − e}), , , , , , , and .

Combining the ranges of V, φ, and [Ca], suitable choices for the voltage, magnetic flux, and calcium scales are and Q = 1 μM. We also see that values of m∞(V), mp∞(V), n∞(V), h∞(V), f([Ca]), G(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 gmax = max{Iextz} = 50 μA. Specifically, Figures 1(a) and 1(b) show the plots of 1/τ(V) and 1/τ(V) over the range V ∈ [−60, 0], which indicates that T ≈ 1.3 ms−1 and T ≈ 0.0025 ms−1. Similarly, we obtain G ≈ 0.06 and G ≈ 1000 pL · ms−1 from Figures 1(c) and 1(d), respectively. So, we have Pmax ≈ 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.

Table 1

Parameter values used in the paper.

Parameter	Value	Parameter	Value	Parameter	Value
C	21 μF	θ n	−29 mV	K CAN	0.74 μM
σ	0.185	θ h	−48 mV	n CAN	0.97
g L	2.8 nS	θ m	−34 mV	[IP3]	1.2 μM
g K	4 nS	θ mp	−40 mV	L IP3	0.37 pL·s−1
g Na	10 nS	σ n	−4 mV	P IP3	31,000 pL·s−1
g NaP	2.8 nS	σ h	6 mV	[Ca]Tot	1.25 μM
g tonic−e	0.3 nS	σ m	−5 mV	K Ca	0.000025 pL−1
g can	0.7 nS	σ mp	−6 mV	V SERCA	400 aMol·s−1
V Na	50 mV	τ¯_n	10 ms	K ERCA	0.2 μM
V K	−85 mV	τ¯_h	10,000 ms	A	0.005 μM·s−1
V L	−65 mV	K I	1.0 μM	K d	0.4 μM
V syn−e	0 mV	K a	0.4 μM

Figure 1

Change of functions of (a) 1/τ(V) with membrane potential V, (b) 1/τ(V) with membrane potential V, (c) G3([Ca]) with [Ca], and (d) gSERCA([Ca]) with [Ca].

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/gmax = 0.42 ms ~ O(1) ms, 1/T = 0.77 ms ~ O(1) ms, 1/T = 400 ms ~ O(100) ms, σ/(Pmax · KCa) ≈ 3.98ms ~ O(10)ms, and 1/Q · A = 200 ms ~ O(100) ms. We use the notation O to denote an order of magnitude estimate: x ~ (10), where n is the nearest integer to log(x). We choose the fast timescale as our reference time, i.e., pick Q = 1 ms, and set

As a result, dimensionless system (7) becomes system (9), namely, with relative rates of all variables: R = O(1), R = O(1), R = O(100), R = O(1), R = O(10), and R = O(100).

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 Ca2+ concentration [Ca] and IP3 channel gating variable l (Equations (3a) and (3b)). For simplicity, we chose gCAN = gCANf([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.

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

Figure 2

The dynamic analysis of system (1a)–(3b) without external stimulations (Iapp = 0). (a) The time course of membrane potential V (black curve). The time course of intracellular calcium concentration [Ca] (red dashed line) is also superposed; (b) one-parameter bifurcation analysis of membrane potential; (c) two-parameter bifurcation analysis of membrane potential.

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 considered to be a constant which we take its average value as h = 0.2834. gCAN 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 gCAN, 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 F1 and F2 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 (gCAN,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 (gCAN,h)-plane. With the increase of gCAN, the trajectory passes through different bifurcation curves. Then, with the decrease of gCAN, 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 Iextz. For other parameter values, please refer to the appendix.

Figure 3

The full system (1a)–(3b) exhibits different firing patterns under different electric currents at k1 = 0.1, k2 = 3 s−1, α = 1 MΩ−1, β = 0.00006 MΩ−1 · V−2 · s−2. The red dashed line represents the time course of intracellular calcium concentration [Ca]. (a) Iextz = −2 μA/cm2; (b) Iextz = 5 μA/cm2; (c) Iextz = 25 μA/cm2; (d) Iextz = 30 μA/cm2; (e) Iextz = 40 μA/cm2; (f) Iextz = 50 μA/cm2.

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.

Figure 4

(a) Partial enlarged area of Figure 3(a). Yellow symbols mark key points along the solution trajectory (star: start of phase ①; triangle: start of phase ②; square: start of phase ③; diamond: start of phase ④); (b) the bifurcation diagram of ISI versus Iextz.

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 increases of current Iextz at first, it begins to decrease rapidly when the current value reaches 21.5 μA/cm2. 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/cm2, the original low potential resting state disappears and all become spiking. In the process of current Iextz increase, phase ③ maintains small amplitude spikings and the duration of time gradually decreases. When the value of Iextz reaches 47.7 μA/cm2, 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/cm2 of an ISI sequence.

3.3.1. Fast-Slow Decomposition

We can obtain one-parameter 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 S-shaped 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 gCAN, 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 F1 and F2 refer to fold bifurcation, and the point HC represents homoclinic bifurcation.

Figure 5

Fast-slow bifurcation analysis of neuronal firing patterns under different electric currents: (a) Iextz = −2 μA/cm2; (b) Iextz = 5 μA/cm2; (c) Iextz = 25 μA/cm2; (d) Iextz = 30 μA/cm2; (e) Iextz = 40 μA/cm2; (f) Iextz = 50 μA/cm2. The black solid line represents the stable nodes (lower branch) and stable focus (upper branch); the black dash-dot line represents the saddles (middle branch) and the unstable focus (upper branch); the red dashed and solid lines represent the maximum and minimum values of the unstable and stable limit cycles, respectively. The points F(i = 1, 2) and subH represent the fold bifurcation and Hopf bifurcation of equilibrium; the points LPC, HC, and SNIC represent the fold bifurcation of the limit cycle, the homoclinic bifurcation, and saddle-node bifurcation on an invariant circle, respectively. The green curve represents the trajectory of the full system.

For Iextz = −2 μA/cm2 and Iextz = 5 μA/cm2, the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable gCAN 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 (F1). 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.

When Iextz increases to 25 μA/cm2 and 30 μA/cm2, the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable gCAN 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 F1. 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 Iextz continuously increases to 40 μA/cm2 and 50 μA/cm2, the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable gCAN 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 gCAN decreases. Take Iextz = 30 μA/cm2 (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 gCAN([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 gCAN([Ca]SNIC) = gCANf([Ca]SNIC).The corresponding periodic orbit in ([Ca], l)-space is shown in Figure 6(b). The point P1 (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 gCAN 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 (P2), 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).

Figure 6

(a) The period of limit cycle of the fast subsystem with respect to parameter gCAN at Iextz = 30 μA/cm2; (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.

3.3.2. Bifurcation Analysis

To better understand the mechanism of the different firing patterns evoked by the electrical current, we present two-parameter bifurcation analysis of the fast subsystem as shown in Figure 7, in which gCAN and h are bifurcation parameters.

Figure 7

Two-parameter bifurcation diagram of the fast subsystem under the electric current in (gCAN,h)-plane showing the curves of fold bifurcation (fc), subcritical Hopf bifurcation (hc), and fold bifurcation of limit cycle (lc), together with projection of the trajectory of the full system (green-colored curve). (a) Iextz = −2 μA/cm2; (b) Iextz = 5 μA/cm2; (c) Iextz = 25 μA/cm2; (d) Iextz = 30 μA/cm2; (e) Iextz = 40 μA/cm2; (f) Iextz = 50 μA/cm2.

For Iextz = −2 μA/cm2 and Iextz = 5 μA/cm2, 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 Iextz increases to 25 μA/cm2 and 30 μA/cm2, 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 Iextz continuously increases to 40 μA/cm2 and 50 μA/cm2, 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 (gCAN, h)-plane changes with the increase of Iextz. That means the value of Iextz affects the critical bifurcations that determine the pattern of bursting. With the increase of Iextz, 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).

3.4. 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 k1. Firing activities of the system corresponding to different values of k1 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 k1 increases.

Figure 8

Firing pattern changes for different parameter values of k1 with Iextz = 0 μA/cm2, k2 = 3 s−1, α = 1 MΩ−1, β = 0.00006 MΩ−1 · V−2 · s−2. The black solid and red dashed curves represent the time course of membrane potential V and intracellular calcium concentration [Ca], respectively. (a) k1 = 0.1; (b) k1 = 0.8; (c) k1 = 1.0; (d) k1 = 1.5.

The ISI bifurcation diagram of the magnetic flow feedback coefficient k1 is shown in Figure 9. Circled numbers are the same as that in Figure 4(a). As k1 increases, the duration of phase ① first slowly increases and then decreases. When k1 increases to 0.88, the original low potential resting state turns to be spiking. With the increase of k1, phase ③ exhibits small amplitude oscillating, and its duration gradually decreases. When the value of k1 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.

Figure 9

The bifurcation diagram of ISI versus k1.

Similarly, Equations (1a), (1b), and (1d) form the fast subsystem, and [Ca] is still the slow variable. We use gCAN = gCANf([Ca]), instead of [Ca], as the bifurcation parameter, and set h to be constant. The bifurcation diagrams of the fast subsystem on the (gCAN,V)-plane are shown in Figure 10. With the parameter k1 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 gCAN decreases. The analysis is similar to that in Figure 6.

Figure 10

Fast-slow bifurcation analysis of neuronal firing patterns with different magnetic flow for (a) k1 = 0.1; (b) k1 = 0.8; (c) k1 = 1.0; (d) k1 = 1.5.

Two-parameter bifurcation of the fast subsystem in (gCAN,h)-plane is shown in Figure 11. Similarly, with the increase of k1, 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 k1 and the electric current have the same effect on the bursting transition.

Figure 11

Two-parameter bifurcation diagram with the effect of magnetic flow on (gCAN,h)-plane for (a) k1 = 0.1; (b) k1 = 0.8; (c) k1 = 1.0; (d) k1 = 1.5.

With the increase of k1, the middle branch of the S-shaped 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 k1 does not change the bifurcation structure but changes the position of the bifurcation point.

Figure 12

The effect of k1 on the bifurcation diagram. The black curve is the bifurcation curve when k1 = 0 and Iextz = 0, that is, the bifurcation curve when there is no external stimulus. The color curves are the bifurcation curves with a nonzero external stimulus. Corresponding values of k1 are marked in the upper right corner of the figure.

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 low-potential resting state to high-potential resting state.

Figure 13

(a) The bifurcation diagram of ISI versus α; (b) the bifurcation diagram of ISI versus β.

4. 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 Iextz 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 Iextz, which indicates that the control of direct current Iextz can lead to desired firing patterns of the system. The variation of magnetic flow can result in a similar effect. The 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.