Abstract
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 -dependent somatic bursting and -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 is the membrane potential and and are gating variables for the voltage-gated potassium and sodium channels, respectively. The functions , , , , , and represent fast sodium current, delayed rectifier potassium current, leakage current, persistent sodium current, calcium-activated nonspecific cationic current, and external tonic drive current, respectively. is a direct current of external stimuli. is the whole cell capacitance. Particularly, The calcium dynamics is given as [13] where represents the fraction of 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 () and the flux out of the cytosol into the ER (). These fluxes are regulated by the intracellular concentration of , [], and channel gating variable, .
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: , where and are fixed parameters [29, 30]. and show the relationship between the membrane potential and the magnetic flux. And we suppose that in Equation (1a). The term 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 lies between 0.95 μM and 1.4 μM [13]. Therefore, in this study, we set μM, which implies that 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 , respectively, such that
Note that , , and are already dimensionless in Equations (1b), (1c), and (3b).
Since the membrane potential typically lies between and , we define over the range and then define , a rescaled version of , by for . As typically lies between 0 μM and 1 μM, we define , over the range . Furthermore, we set and and then define . We also define
According to system (1a)–(3b), we get the following dimensionless system: with dimensionless parameters , , , , , , , , and .
Combining the ranges of , , and , suitable choices for the voltage, magnetic flux, and calcium scales are and . We also see that values of , , , , , , , , , and all lie in the range . Combining the parameter values in Table 1 and the values of variables in the paper, we have . Specifically, Figures 1(a) and 1(b) show the plots of and over the range , which indicates that and . Similarly, we obtain and from Figures 1(c) and 1(d), respectively. So, we have . Using these values, we see that all terms in the right-hand sides of Equations (7a)–(7f) are bounded (in absolute value) by one.
(a)
(b)
(c)
(d)
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 , , , , and . We use the notation to denote an order of magnitude estimate: , where is the nearest integer to . We choose the fast timescale as our reference time, i.e., pick ms, and set
As a result, dimensionless system (7) becomes system (9), namely, with relative rates of all variables: , , , , , and .
From this, we conclude that , , and evolve on a fast timescale, evolves on a slow time scale, and and 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 channel gating variable (Equations (3a) and (3b)). For simplicity, we chose as the bifurcation parameter, where is a monotonically increasing concave function of [Ca].
According to the nondimensionalization of the model, the variable is super slow. So, can be considered a constant when we do the fast-slow decomposition. That is, we take 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 (black solid) and the intracellular calcium concentration [Ca] (red dashed) are shown in Figure 2(a).
(a)
(b)
(c)
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 is a super slow variable. So, is considered to be a constant which we take its average value as . 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 , 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 and 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 (,)-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 (,)-plane. With the increase of , the trajectory passes through different bifurcation curves. Then, with the decrease of , 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 . 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.
(a)
(b)
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 at first, it begins to decrease rapidly when the current value reaches 21.5 . 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 , the original low potential resting state disappears and all become spiking. In the process of current increase, phase ③ maintains small amplitude spikings and the duration of time gradually decreases. When the value of reaches 47.7 , 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 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 , 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 and refer to fold bifurcation, and the point HC represents homoclinic bifurcation.
(a)
(b)
(c)
(d)
(e)
(f)
For and , the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable is shown in Figures 5(a) and 5(b), respectively. Here, we set and , respectively ( takes the average value). Firstly, the trajectory transits from the lower rest state to the spiking state via fold bifurcation (). 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 increases to 25 and 30, the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable is shown in Figures 5(c) and 5(d), respectively. In these two cases, the parameter is set to be and , respectively. The trajectory jumps up to the upper branch of the S-shaped curve via fold bifurcation . 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 continuously increases to 40 and 50 , the bifurcation of fast subsystem (1a), (1b), and (1d) with respect to the slow variable is shown in Figures 5(e) and 5(f), respectively. is set to be and , 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 decreases. Take (Figure 5(d)) as an example, the periodic orbit of the (,)-space is shown in Figures 6(a) and 6(b). Let be the value of [Ca] at SNIC; the period of the limit cycle tends to infinite when at which the saddle-node bifurcation on an invariant circle (SNIC) occurs. Then, we get according to .The corresponding periodic orbit in (, )-space is shown in Figure 6(b). The point (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 correspondingly, and this increase pushes the trajectory into the branch of periodic orbits. Then, the value of [Ca] slowly decreases while the trajectory in (, )-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 (), the system keeps firing. Once [Ca] falls below , the trajectory solution stops firing and jumps down to the lower branch of the S-shaped curve (Figure 5(d)). Next, the trajectory in (, )-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).
(a)
(b)
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 and are bifurcation parameters.
(a)
(b)
(c)
(d)
(e)
(f)
For and , as shown in Figures 7(a) and 7(b), the trajectory of the full system passes through the fold bifurcation curve (), the Hopf bifurcation curve (), and the fold limit cycle bifurcation curve () at the same time. The fold bifurcation curve and the Hopf bifurcation curve intersect near the full system trajectory.
When increases to 25 and 30, the trajectory of the full system moves down with the decreasing of . It still passes through the three bifurcation curves of fc, hc, and lc. But the gap distance between the fold bifurcation () and the Hopf bifurcation curve () increases, as shown in Figures 7(c) and 7(d).
When the value of continuously increases to 40 and 50 , the trajectory of the full system continues to move down to the decreasing direction of . It passes through two bifurcation curves of and . 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 (, )-plane changes with the increase of . That means the value of affects the critical bifurcations that determine the pattern of bursting. With the increase of , 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 . Firing activities of the system corresponding to different values of 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 increases.
The ISI bifurcation diagram of the magnetic flow feedback coefficient is shown in Figure 9. Circled numbers are the same as that in Figure 4(a). As increases, the duration of phase ① first slowly increases and then decreases. When increases to 0.88, the original low potential resting state turns to be spiking. With the increase of , phase ③ exhibits small amplitude oscillating, and its duration gradually decreases. When the value of 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 , instead of [Ca], as the bifurcation parameter, and set to be constant. The bifurcation diagrams of the fast subsystem on the (,)-plane are shown in Figure 10. With the parameter 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 decreases. The analysis is similar to that in Figure 6.
(a)
(b)
(c)
(d)
Two-parameter bifurcation of the fast subsystem in (,)-plane is shown in Figure 11. Similarly, with the increase of , 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 and the electric current have the same effect on the bursting transition.
(a)
(b)
(c)
(d)
With the increase of , 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 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 low-potential resting state to high-potential resting state.
(a)
(b)
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 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 , which indicates that the control of direct current 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.
Appendix
For , the function of takes the form .
For , the function of takes the form .
The activation of the CAN current by the calcium concentration is given as
[], as well as , is described as follows: in which is given as
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest in this work.
Acknowledgments
The authors thank Dr. Yangyang Wang for valuable suggestions. This work is supported by the National Natural Science Foundation of China (Grant No. 11872003) and North China University of Technology Research Fund Program for Key Discipline (No. 110052972027/014).