Pre-Bötzinger complex is considered to have a closely relationship with the respiratory rhythms. In this paper, we investigate the bursting phenomena of the pre-Bötzinger complex respiratory neuron with periodic slow stimulation. Since the external forcing visit the spiking and rest areas in different ways, thus the system could generate various complex bursting patterns. With the external forcing is taken as a slow variable that modulates the dynamics of the system, different types of bursting are distinguished and the generation mechanism is explored by using the combination of two parameter bifurcation analysis and fast slow decomposition. Our results show that both the external forcing and the neural intrinsic property play an important role in neural activities.

## I. INTRODUCTION

Respiration is an important sign of life. It is known that the pre-Bötzinger complex in mammalian brain stem plays a crucial role in the generation of rhythmical respiratory.^{1} Within pre-Bötzinger complex, there exist many neurons which can yield different firing activities, including quiescence, tonic spiking and bursting.^{2} Here we are mainly concerned with bursting activities, which have a closely relationship with the generation of respiratory rhythm.^{3} Due to the importance of bursting, much work has been done to investigate the bursting generation mechanism,^{4–6} and fast-slow analysis is a powerful tool in the study of bursting activities.^{7} But, the following two aspects had not been fully considered in the previous studies. On the one hand, the influence of external stimulus is rarely considered, however, such effects are factual.^{8} On the other hand, most previous studies had largely focused on the bursting phenomena of autonomous systems, but the work on nonautonomous system has been rarely reported.^{9,10} Taking these processes into account, we consider a pre-Bötzinger complex respiratory neuron model with periodic external stimulation, and investigate the influence of persistent sodium current on the behavior dynamics.^{4} Appropriate parameters are taken such that there exists an order gap between the natural frequency and the forcing frequency,^{9,10} therefore the model is characterized by the multi-time scale feature, which provides the possibility of the appearance of the bursting behavior. The parameter regions for different bursting patterns in the model driven by the external forcing are presented and classified by using both two parameter bifurcation and fast slow dynamics analysis.

This paper is organized as follows. In Section II, we describe the pre-Bötzinger respiratory neuron model with external excitation. Two parameter bifurcation analysis is given in Section III, and bursting transitions from one bursting pattern to another are presented numerically in Section IV. In the calculation process, the fourth order Runge-Kutta method with the step length 0.01 is used, and the numerical simulations have been achieved by XPP^{11} and MATLAB. The conclusion is given in the last section.

## II. MODEL DESCRIPTION

In order to discuss the complex bursting pattern for the higher dimensional system, a periodic stimulus term is added into the pre-Bötzinger complex respiratory neuron proposed by Butera *et al* in 1999,^{11} thus a new three-dimensional nonautonomous system is produced, where *A* and *w* represent the amplitude and frequency of the external forcing, respectively. This modified model is described as follows:

Where the first equation is the dynamics about the membrane potential *v*, *C* is the cell capacitance, and *t* is the time. *I*_{Nap}, *I*_{Na}, *I*_{K} and *I*_{L} represent persistent sodium current, fast sodium current, delayed-rectifier potassium current and leak current, respectively. The expressions for these currents are modeled by the Hodgkin-Huxley formulation, more specifically:

The second and the third equations are the dynamics for the gating variables *h* and *n*, respectively, where *x*_{∞}(*v*) ($x\u2208{mp,mNa,h,n}$)is the steady-state voltage-dependent (in)activation function and *τ*_{x}(*v*) is the voltage-dependent time constant, they are described according to

Here the external forcing frequency w is chosen as 0.0001*Hz*, thus this very small number enables the term *A*^{*}sin(*ωt*) to vary slowly, thereby it can be considered as a slow-varying parameter for the system. As a result, we can perform the fast-slow analysis to study the effect of the slow forcing on the dynamics of pre-Bötzinger neuron, where *A*^{*}sin(*ωt*) can be seen as a slow variable, while *V*, *n* and *h* are fast variables. In the following study, we choose the persistent sodium current conductance *g*_{Nap} as a control parameter, other parameter values are fixed as *C* = 21 *pF*, *g*_{Na} = 1.5*nS*, *g*_{K} = 11.2*nS*, *g*_{L} = 3*nS*, *E*_{Na} = 50*mV*, *E*_{K} = −85*mV*, *E*_{L} = −55*mV*, *τ*_{n} = 10*ms*, *τ*_{h} = 100*ms*, *σ*_{p} = −6, *σ*_{m} = −5, *σ*_{h} = 6, *σ*_{n} = −3, *θ*_{p} = −40, *θ*_{m} = −34, $\theta h=\u221248$, $\theta n=\u221225$.

## III. TWO PARAMETER BIFURCATION ANALYSIS

Two parameter bifurcation analysis of the whole system with respect to the slow forcing *A*^{*}sin(*ωt*) and control parameter *g*_{Nap} is illustrated in Fig. 1 (a). For clarity, here only the pertinent details are presented. Fig. 1 (b) is the local magnification of the rectangle in Fig. 1(a). The codimension-1 bifurcation curves in Fig. 1 (a) are composed of two fold bifurcation curves f1and f2 (blue lines), two Hopf bifurcation curves h1 and h2 (black lines), two homoclinic bifurcation curve hc1 and hc2 (pink lines), and a fold limit cycle bifurcation curve lc (red line). Besides these codimension-1 bifurcations, the system also undergoes three codimension-2 bifurcation points, that is, two Bogdanov-Taken bifurcation points (BT_{1} and BT_{2}) and a Cusp bifurcation point (CP), which can be seen clearly in Fig. 1(b). Here the amplitude of the periodic stimulation is taken as *A* = 50 such that all the codimension-1 curves could be involved in the oscillation area of the slow forcing *A*^{*}sin(*ωt*), see Fig. 1(a).

We will explain that both the stable limit cycle and Hopf bifurcations play a crucial role in deciding the bursting type. To explain more fully, it is desirable to use A to denote the intersection of three curves f1, lc and hc2, where the bifurcation curves in the region nearby can be seen more clearly in Fig. 1(b). In addition, two inflection points of the Hopf curve h1 are represented by B and C, respectively. In fact, point A gives the stable limit cycle an opportunity to participate in the bifurcation of the system, which is directly related to the generation of bursting. Inflection points B and C determine the number of Hopf bifurcations, which could lead the system to exhibit different bursting rhythms. Thus three points A, B and C can be seen as the dividing points to separate different bursting patterns. By taking some representative values, four typical bursting oscillations are illustrated in Fig. 2.

## IV. THE GENERATION MECHANISM OF DIFFERENT BURSTING TYPES

In this section, we will investigate the generation mechanism of complex bursting patterns illustrated in Fig. 2. In order to better reveal the dynamic behavior, time evolution of the slow forcing 50^{*}sin(0.0001*t*) are plotted in Figs. 3(a)–6(a), and the corresponding bifurcation diagram are given in Figs. 3(b)–6(b).

### A. The case of *g*_{Nap} = 1

In this case, system exhibits the bursting pattern as shown in Fig. 2(a). From Fig. 1 (b), it can be seen that *g*_{Nap} = 1 is located below the critical point A, which indicates that except limit cycle bifurcation, all other bifurcations including fold bifurcations F1 and F2, hopf bifurcations H1 and H2, as well as homoclinic bifurcations HC1 and HC2 are involved in the bifurcation diagram, see Fig. 3(b). It is found that the equilibrium points of the system with respect to the slow forcing forms a *S*-shaped curve, which is made up of three branches. By computing the eigenvalues of the Jacobian matrices at the equilibriums, it is concluded that the lower branch is consisted of saddle-focus which possess property both of the node and the focus together. The middle branch is composed of unstable saddles, and the upper branch is made up of focus. The stable limit cycles originated from Hopf bifurcation H1 (H2) terminate at the homoclinic bifurcation HC1 (HC2) in the middle branch. The phase portrait of the system is also imposed in Fig. 3 (b).

The dynamics of the system can be explained by use of the *S*-shaped curve and the phase portrait together. In order to better understand the bursting pattern, we start to explain the system behavior at the silent phase. It is seen that the trajectory firstly follows the lower branch until the stable state become unstable via Hopf bifurcation H2, which means that the system should come into the spiking state. But due to the slow passage effect^{12} the trajectory does not jump up immediately, instead, the trajectory still lingers on the lower branch until the fold bifurcation F2 is encountered. Note that the slow forcing is a periodic function, thus its value begins to decrease after it reaches its maximum. Correspondingly the trajectory switches to left and jump down via the fold bifurcation F1 rather than Hopf bifurcation H1, and slow passage effect is also responsible for this fact. Finally the phase portrait returns back to begin another period. Both the disappearance of the silent phase and active phase depend on the Hopf bifurcations. In addition, as described above, there exists a hysteresis loop of “fold/fold” type. Two fold bifurcations F1 and F2 correspond to the turning points of the firing state, which are displayed with two red lines in Fig. 3(a). The region between the red lines is spiking area (SA) of bursting, and the remaining regions are resting areas. Specially, the resting state corresponding to the lower branch is denoted by RA1, and the resting state corresponding to the upper branch is denoted by RA2. According to the classification of bursting type,^{7} this type of bursting is referred to “Hopf/Hopf” bursting via “fold/fold” hysteresis loop.

### B. The case of *g*_{Nap} = 2

From Fig. 1(b) it can be seen that now the limit cycle bifurcation is involved in this case, which is crucial for complex bursting generation. The periodic movement with mixed mode oscillations is given in Fig. 2(b). From the bifurcation diagram Fig. 4(b), it is seen that the unstable limit cycle originated from H1 become stable via the fold limit cycle bifurcation LC, and finally terminates at Hopf bifurcation H2. In succession, the dynamic behavior of bursting is discussed. The trajectory firstly follows the lower branch which corresponds to the resting state of bursting, and then switches to the stable limit cycle, thus the system gets into the spiking state until the track experiences the second Hopf bifurcation H2. The trajectory keeps on move rightward until the sine excitation reaches its maximum. After that, the portrait turns to the left until returns to Hopf bifurcation H2, and the subsequent behavior is similar to what had happened on the right side of H1 described above. As the former case, two Hopf bifurcations are responsible to the beginning and the end of the bursting, thus this type of bursting is also defined as “Hopf/Hopf” bursting.^{7} The oscillation region of the slow variation are divided into three parts (two resting areas RA1 and RA2, and a spiking area SA) by two Hopf bifurcations rather than two fold bifurcations, which is the major difference from the previous situation. Furthermore, since the slow variable visits the spiking area SA twice time in one period, and thus the bursting is composed of two sub-burstings.

### C. The case of *g*_{Nap} = 3

In this case, the system generates a new type of bursting, as shown in Fig. 2(c). From Fig. 1 (a), the parameter value is located in the region between two inflection points B and C. Quite evidently, now the bifurcation curve has four Hopf bifurcations, see Fig. 5(b), thus the oscillating region of the slow forcing is divided into three resting areas (RA1, RA2 and RA3) as well as two spiking areas (SA1 and SA2), see Fig. 5(a). Similar fast slow analysis is used to analyze this type of bursting. The bursting trajectory is located at the lower branch firstly, which means that slow variable is in the resting area RA1. Then it passes through the first Hopf bifurcation H1 into the spiking area SA1, and starts to wrap towards the stable limit cycle until the second Hopf bifurcation H2 is encountered, leading the slow forcing introduced in the resting area RA2. As time progresses, the trajectory undergoes the third Hopf bifurcation H3, resulting in the appearance of the stable limit cycle, which indicates that the slow variable enters into the spiking area SA2, thus the system is in the process of firing state. Finally, the trajectory passes through the fourth Hopf bifurcation H4, and correspondingly the slow variable would enter into the resting area RA3. After a while, the slow forcing oscillates to the left and goes through the reverse order, that is, resting area RA3-spiking area SA2-resting area RA2-spiking area SA1-resting area RA1, thereby a bursting cycle is completed. It is seen that this type of bursting still belongs to “Hopf/Hopf” type of bursting.^{7} The key difference is that the during a period the slow forcing visits the spiking areas SA1 and SA2 on four occasions, thus the bursting is composed of four sub-burstings.

### D. The case of *g*_{Nap} = 4

From two parameter bifurcation diagram Fig. 1(a), it can be seen that now the parameter value is just located above the convex point C, which indicates that the system has only experienced two Hopf bifurcations, and fast slow dynamics analysis Fig. 6(b) makes this fact evidently. The generation mechanism of this type of bursting is similar to that of the case with *g*_{Nap} = 1. It is seen that trajectory is located in the silent phase of bursting firstly. And then the trajectory is switched to the active phase when Hopf bifurcation H1 is reached. The trajectory starts wrap around the stable limit cycle and terminates at the Hopf bifurcation H2, which means that the repetitive spiking state is about to disappear. In what follows, the slow forcing begins to decrease after it reaches the maximum. Thus the trajectory turns to the left until Hopf bifurcation H2 is encountered. The behavior of the orbit on the left side of H2 is similar to that of the orbit on the right side of H1described above. According to the classification of the bursting type,^{7} this type of bursting is classified as “Hopf/Hopf” type of bursting. The slow forcing visits the spiking area SA twice time in a period, thus the bursting is composed of two sub-burstings. In addition, Compared with Fig. 3(a), it is seen that in this case the two Hopf bifurcations are positioned at a greater distance, which leads to the long duration in the repetitive spiking, see Fig. 2(d).

## V. CONCLUSION

In summary, pre-Bötzinger respiratory neuron with slow periodic forcing could exhibit complex bursting patterns with different waveforms. The mechanisms underlying the appearance of these burstings can be understood by fast slow dynamics analysis. In fact, the variation region of the slow forcing can be divided into resting areas and spiking areas by bifurcation points. Bursting oscillations are created since the slow forcing visits the resting areas and the spiking areas alternately. The bursting patterns can have different clusters of repetitive spiking, such as two sub-burstings and four sub-burstings, which are decided by the visiting modes that the slow forcing exhibits during its evolution. The boundaries of different bursting patterns are obtained by two parameter bifurcation analysis. Here we point out that if we do not consider the outside periodic control, the original model will be in the quiescent or infinite long terms with no-firing state, which seems could be called as “respiratory death”, thus the external forcing is necessary for bursting generation. Therefore, this study shows that the periodic forcing has great influence on the neuron firing activities, and provides basis for further research on the generation of the respiratory rhythm.

## ACKNOWLEDGMENTS

This work is supported by the National Natural Foundation of China (Nos. 11402057 and 11872084) and Youth Innovation Talents Project of Colleges and Universities in Guangdong Province (No. 2014KQNCX137), Guangdong province Medical Science and Technology Research Fund Project (No. A2016147), Guangdong Province Data Science and Engineering Research Center’s Open Fund Project (No. 2016KF02), Natural Science Foundation for the Higher Education Institutions of Anhui Province of China (Nos. KJ2016SD54 and KJ2017A460).