Fast–slow variable dissection with two slow variables related to calcium concentrations: a case study to bursting in a neural pacemaker model

Neuronal bursting is an electrophysiological behavior participating in physiological or pathological functions and a complex nonlinear behavior alternating between burst and quiescent state modulated by slow variables. Identification of dynamics of bursting modulated by two slow variables is still an open problem. In the present paper, a novel fast–slow variable dissection method with two slow variables is proposed to analyze the complex bursting simulated in a four-dimensional neuronal model to describe firing associated with pathological pain. The lumenal (Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document}) and intracellular (Cin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{in}$$\end{document}) calcium concentrations are the slowest variables, respectively, in the quiescent state and burst duration. Questions encountered when the traditional method with one slow variable is used. When Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document} is taken as slow variable, the burst is successfully identified to terminate near the saddle-homoclinic bifurcation point of the fast subsystem and begin not from the saddle-node bifurcation. With Cin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{in}$$\end{document} chosen as slow variable, Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document} value of the initiation point of burst is far from the saddle-node bifurcation point, due to Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document} not contained in the equation of the membrane potential. To overcome this problem, both Cin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{in}$$\end{document} and Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document} are regarded as slow variables; the two-dimensional fast subsystem exhibits a saddle-node bifurcation point, which is extended to a saddle-node bifurcation curve by introducing the Clum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_\mathrm{lum}$$\end{document} dimension. Then, the initial point of burst is successfully identified to be near the saddle-node bifurcation curve. The results present a feasible method for fast–slow variable dissection and a deep understanding of the complex bursting behavior with two slow variables, which is helpful for the modulation to pathological pain.

est variables, respectively, in the quiescent state and burst duration. Questions encountered when the traditional method with one slow variable is used. When C lum is taken as slow variable, the burst is successfully identified to terminate near the saddle-homoclinic bifurcation point of the fast subsystem and begin not from the saddle-node bifurcation. With C in chosen as slow variable, C lum value of the initiation point of burst is far from the saddle-node bifurcation point, due to C lum not contained in the equation of the membrane potential. To overcome this problem, both C in and C lum are regarded as slow variables; the two-dimensional fast subsystem exhibits a saddle-node bifurcation point, which is extended to a saddle-node bifurcation curve by introducing the C lum dimension. Then, the initial point of burst is successfully identified to be near the saddlenode bifurcation curve. The results present a feasible method for fast-slow variable dissection and a deep understanding of the complex bursting behavior with two slow variables, which is helpful for the modulation to pathological pain.

Introduction
Identification of nonlinear dynamics of neuronal electrical behavior is very important for understanding the functions of nervous system [1][2][3][4][5][6]. Neuronal bursting is a complex nonlinear electrical behavior related to fast and slow variables, which participates many kinds of physiological functions, pathological functions, and even brain diseases [7][8][9]. For example, bifurcations of bursting have been identified to participate information coding [10][11][12][13][14] in blood pressure and pathological pain. In an animal model of injured nerve fiber to investigate the pathological pain [15], or called neuronal pacemaker, many complex bursting patterns containing periodic, chaotic, and stochastic patterns have been observed in the electrophysiological experiment [16]. Modulated by potassium or calcium concentrations, the bursting patterns exhibit complex bifurcations, which is related to information coding of firing patterns instead of firing frequency [12,13]. Recently, different bursting patterns have built relationships to the spreading depression (SD) and seizure (SZ) [17][18][19][20] and respiratory function [21][22][23][24]. The bursting in the lateral habenula (LHb) has been related to depression [25]. Bursting activity related to autapse in the pyramidal neuron is associated with coincidence detection [26]. Identification of bursting dynamics has been an important issue of nonlinear dynamics [27][28][29][30].
Typical bursting exhibits alternation behavior between burst and quiescent state, which is modulated by fast and slow variables [1,22,[31][32][33][34][35]. In neurophysiology, the fast variables are mainly related to the membrane potential and sodium or potassium channel, and the slow variables are often associated with calcium concentrations. The minimal theoretical model to describe bursting is always composed of three variables. For example, in the famous Chay model and Hindmarsh-Rose model, except for a slow variable, there are two fast variables comprising the fast subsystem whose differential equations contain the slow variable [36][37][38]. The bursting dynamics can be well explained with a feasible method which is called the fast-slow variable dissection [1,39]. With the slow variable regarded as bifurcation parameter, the fast subsystem often exhibits spiking behavior and resting state and bifurcations between the two behaviors. The bursting behavior is determined by the dynamical behaviors and bifurcations of the fast subsystem. In general, the burst corresponds to the spiking behavior and the quiescent state to the resting state. Furthermore, the transitions between burst and quiescent state of bursting are associated with the bifurcations between the two behaviors of the fast subsystem, spiking and resting state. Correspondingly, a bursting pattern is named with the two bifurcations, and many bursting behaviors are classified into different patterns with such a classification method. For example, for the well-known "Fold/Homoclinic" bursting pattern, the former and latter bifurcations are associated with the beginning and ending phases of the burst. The dynamics of bursting with one slow variable has been successfully analyzed with the fast-slow variable dissection [7,[39][40][41][42][43][44][45][46][47].
In many real nervous systems, the bursting patterns are related to two or more slow variables. For example, both intracellular and lumenal calcium concentrations are slow variables for many nervous system such as the injured nerve fiber [14], pyramidal neuron of CA1 region of hippocampus [48], neuron of pre-Bötzinger (PB) complex [21,23,41,[49][50][51][52][53], and β-cell [54]. At least three cases of two slow variable can be found from the theoretical models. For case-1, two slow variables exhibit very different time scales. One is much slower than the other. The representative example refers to the β-cells with two slow variables. Therefore, when the fast-slow variable dissection is performed, the slowest one is taken as a constant value and the other slow one is regarded as the bifurcation parameter. The dynamics of such bursting can be explained to a certain extent. Case-2 corresponds to two variables with approximate time scale. The pyramidal neuron of CA1 region of hippocampus related to SD or SZ [55][56][57] belongs to such case. Both potassium and sodium concentrations are slow variables with nearly approximate time scale. For such a theoretical model, question is encountered when the fast-slow variable dissection with one slow variable is performed. If only one slow variable is taken as bifurcation parameter, the bifurcation points between the spiking and resting state of the fast subsystem mismatch the transition points between the burst and quiescent state of the bursting, as reported in Ref. [58]. To overcome this problem, fast-slow variable dissection method with two slow variables is proposed, with both low variables (potassium and sodium concentrations) regarded as the bifurcation parameters of the fast subsystem. The codimension-1 bifurcation curves in the parameter plane of potassium and sodium concentrations are acquired. The dynamics of bursting is identified according to the codimension-1 bifurcation curves intersected with the phase trajectory of the bursting. The bursting for SD is related to the Hopf bifurcation and for SZ to the saddle-node bifurcation. For the case-2, the equations of the fast subsystem contain both slow variables of potassium and sodium concentrations.
Except for case-1 and case-2, the two slow variables are intracellular calcium concentration (C in ) and lumenal calcium concentration (C lum ) for case-3, and the representative examples contain the model for neuron in pre-Bötzinger (PB) complex [51][52][53] and a four-dimensional theoretical model to describe the firing behavior of the injured nerve fiber. For the fourdimensional theoretical model, the variables V and n are considered to describe the membrane potential and the gating dynamics of potassium channel, respectively, except for the variables C in and C lum . The differential equations of V and n do not contain the variable C lum , and the differential equation of C lum does not contain the variables V or n [51][52][53]. Therefore, the fast-slow variable dissection method with two slow variables for case-1 [54] and case-2 [58] cannot be directly adopted to identify the dynamics of bursting for case-3, due to that the fast subsystem of V and n does not contain the variable C lum . When the fast-slow variable dissection is performed to analyze the dynamics of case-3 bursting with two slow variables, question arises in two aspects.
(1) If C lum is regarded as slow variable, the fast subsystem is composed of the differential equations of V , n, and C in . The differential equations of V and n (not contain the variable C lum ) and the differential equation of C lum (not contain V and n) become uncoupling. How to acquire the dynamics including the bifurcations of the fast subsystem becomes a problem to be answered. (2) If both C lum and C in are taken as slow variables, how to acquire the bifurcations of the fast subsystem composed of differential equations of V and n in two-parameter plane (C in , C lum ) is a problem, due to that the equations of fast subsystem do not contain C lum . Therefore, identification of dynamics of bursting modulated by C in and C lum is still an open question.
To answer the questions mentioned above, the bursting modulated by C in and C lum in the four-dimensional pacemaker model to describe dynamics of the injured nerve fiber is studied in the present paper. The variable C in is identified to be the slowest variable in the quiescent state, and C lum is the slowest variable in the burst duration. A feasible process to perform the fastslow variable dissection with two slow variables C in and C lum is proposed. The dynamics of the bursting is acquired by three manners of the fast-slow variable dissection, with C in , C lum , and both C in and C lum taken as slow variables, respectively. Firstly, when C lum is taken as the slow variable, the termination phase of the burst is identified to correspond to the saddle-homoclinic bifurcation of the fast subsystem (V , n, C in ), while the initial phase mismatches the saddle-node bifurcation, due to that C lum instead of C in is the slowest variable in the burst duration. Secondly, when C in is chosen as the slow variable, the initiation phase of burst is close to the saddle-node bifurcation point of the fast subsystem (V , n, C lum ) in C in direction but far from the saddle-node bifurcation point in C lum direction, due to the uncoupling between the differential equations of V and n and the equation of C lum . Thirdly, to overcome these problems, both C in and C lum are regarded as the slow variables. Due to that the saddle-node bifurcation with respect to C in is independent of C lum (the equations of V and n not contain C lum ) in the fast subsystem, the saddle-node bifurcation point is extended to a saddle-node bifurcation curve by introducing the dimension of C lum , with the values of the saddle-node bifurcation point along the C lum -axe unchanged. Then, the burst is identified to begin from a point near the saddle-node bifurcation curve, and the C lum value of the begging point of the burst is determined by not the fast subsystem but the full system. Last, a comprehensive viewpoint of dynamics of the bursting is presented by comparisons between dynamics acquired with the three manners of the fast-slow variable dissection. These results present deep understandings of the complex dynamics underlying the bursting behaviors of the neural pacemaker model, which is helpful to modulate the bursting behaviors related to the pathological pain, and provide a novel feasible process to perform the fast-slow variable dissection to bursting patterns, which can be extensively utilized to analyze the bursting of other models with slow variables C in and C lum .
The rest of the present paper is organized as follows. Sections 2, 3, and 4 present the model and method, results, and conclusion, respectively.

Model
The neuronal pacemaker model is proposed based on the results of electrophysiological experiments on the injured nerve fibers (or called neural pacemaker) [11,12,[14][15][16] and the Chay model, which exhibits complex bursting behaviors and bifurcations from bursting to spiking patterns. The differential equations of the four-dimension pacemaker model are described as follows: where the variable V describes the membrane potential, n describes the activation probability of voltageactive potassium K + channel, C in is the intracellular concentration of Ca 2+ , and C lum represents the lumenal calcium concentration in the endoplasmic reticulum. The parameters g i , g kv , g l , and g kc are the corresponding maximal conductance for Na + , voltage-active K + , leakage, and intracellular concentration-dependent K + ions, respectively. And the parameters v i , v k , v l , and v c are the corresponding equilibrium potentials. The functions h ∞ , m ∞ , and n ∞ are the inactivation probability of Na + channel, the activation probability of Na + channel, and the steady-state value of the second variable n, respectively. More details about the model refer to Refs. [14,30,36]. The equations of these functions are given as follows: The unit of membrane potential is mV, the unit of time is ms, and the unit of conductance is mS. In the present paper, v i = 100 mV, v k = −75 mV, v l = −40 mV, g i = 1800 mS, g kv = 1700 mS, g kc = 16 mS, g l = 7 mS, k c = 3.3/18, τ c = 1.0, λ n = 237.3, k erp = 15, k rel = 0.3, and k crc = 0.5. The parameter v c is regarded as the bifurcation parameter.

Methods
The four-dimensional pacemaker model is solved using fourth-order Runge-Kutta integration method, and the integration time step is 0.001 ms. If the upstroke of the membrane potential reaches − 30.0 mV, a spike is thought to appear. In fact, the spike threshold can be chosen as a value between −25 and −42 mV. The bifurcations are acquired using the software XPPAUT 8.071 [59].

Bursting patterns and bifurcations
The theoretical model exhibits various bursting patterns and spiking patterns at different v c values. For example, period-4 bursting appears when v c = 110 mV, as shown in Fig. 1a. The burst and a long quiescent state alternate. The inter-spike intervals (ISIs) of the firing pattern can be defined as follows: where t k is the time of the peak value of the kth spike and N is the number of spikes. In the present paper, N = 200 is used. For the period-4 bursting shown in Fig. 1a, values exhibit 4 fixed values in turn. With decreasing v c , the ISIs of the firing exhibit a complex bifurcation process, as depicted in Fig. 1b. Such a bifurcation process has been observed in the experiments [11,12,14]. The bursting behaviors include the chaotic and periodic bursting patterns such as period-1, -2, . . ., -11, . . ., and even -17. In the present paper, the period-11 bursting for v c = 80 mV is chosen as representative to show the results, as shown in Fig. 2a1. 3.2 C in and C lum are slow variables in quiescent state and burst durations, respectively The variables V , n, C in , and C lum for the period-11 bursting are depicted in Fig. 2a1, b1, c1, and d1, respectively. And the derivatives dV /dt, dn/dt, dC in /dt, and dC lum /dt are shown in Fig. 2a2, b2, c2, d2, respectively. The blue and red stars represent the termination and initiation phases of the burst, respectively. During the quiescent state, dC in /dt (green) is closer to zero than those of dV /dt (black), dn/dt (red), and dC lum /dt (blue), as shown in Fig. 3a. In the burst duration, the value of dC lum /dt (blue) is closer to zero than those of dC in /dt (green), dV /dt (black), and dn/dt (red), as illustrated in Fig. 3b. C in exhibits the slowest variation during the quiescent state, and C lum manifests the slowest variation in the burst duration. Therefore, C in and C lum are, respectively, slow variables during quiescent state and burst within a period of the period-11 bursting.

The question encountered when fast-slow variable dissection is performed
Due to that C in is the slowest variable in the quiescent state and C lum is the slowest variable in burst duration, the fast-slow variable dissection is performed in 3 manners in the present paper: (1) C lum is taken as the slow variable, and equations of V , n, and C in form the fast subsystem; (2) C in is chosen as the slow variable, and equations of variables V , n, and C lum are composed of the fast subsystem; (3) C lum and C in are regarded as the slow variables, and equations of V and n are the fast subsystem.
It should be noticed that a characteristic of the model makes it difficult to perform the fast-slow variable dissection to the bursting behavior. Both of Eqs. (1) and (2) (V and n) do not contain the variable C lum . Both of variables V and n are not contained in the differential equation of C lum (Eq. (4)), and Eq. (4) is only determined by the variables C in and C lum . Therefore, there is no direct relationship between the slow variable C lum and the fast variables V and n. Such a characteristic of the model makes it difficult to perform the fast-slow variable dissection for the manners 2 and 3, which is introduced when the problem is encountered.

Burst terminates at the saddle-homoclinic bifurcation: manner 1
The equations of three variables V , n, and C in (Eqs. (1)-(3)) form the fast subsystem, where C lum is taken as the bifurcation parameter. With respect to C lum , the equilibrium point curve of the fast subsystem (V , n, C in ) is called curve E1 in the present paper, as depicted in Fig. 4. The projection of E1 in (C lum , C in , V ) and (C lum , C in , n) space is illustrated in Fig. 4a, b, respectively. The projection of E1 exhibits a Z-shaped curve containing three branches. As C lum decreases, the middle branch which represents a saddle and the lower one representing a stable node intersect with each other to form a saddle-node bifurcation point labeled as SN1 (red diamond). After the bifurcation, both the saddle and the node disappear. For the fast subsystem (V , n, C in ), there exists a stable limit cycle with maximal (LCU1, upper pink curve) and minimal (LCL1, lower pink curve) values, as shown in Fig. 5a. The projection of the equilibrium point curve E1 (red) and the trajectory (cyan) of period-11 bursting on (C lum , C in , V ) space are depicted in  Fig. 5b. When the limit cycle contacts with the middle branch which corresponds to a saddle, a saddle-homoclinic (SH) bifurcation point (blue diamond almost overlapping with the blue star, labeled as SH1) emerges. At the SH1 point, the stable limit cycle becomes a homoclinic orbit begging from and ending at the saddle appears. The stable limit cycle disappears after the SH1 bifurcation. The bursting trajectory (cyan) evolves in clockwise. Then, the burst is identified to terminate at the SH1 point. After the SH1 point, the behavior of bursting runs down to the lower branch of the curve E1, i.e., the stable node, and then runs along the lower branch. After running across the SN1 point, the trajectory goes through a relative long distance, and then runs up to the limit cycle (pink) to form the burst with 11 spikes. Therefore, the burst begins not at the SN1 point but a point (red star) far away from the SN1 point in C lum distance.
Four panels of Fig. 6 show the dynamics depicted in Fig. 5 from different angles of view. The projection of E1 (red), maximal (LCU1, upper pink curve) and minimal (LCL1, lower pink curve) values of the cycle on the (C lum , V ) plane is depicted in Fig. 6a. The projections of the part around the trajectory (cyan) of period-11 bursting on the (C lum , V ), (C in , n), and (C lum , C in ) planes are enlarged, as depicted in Fig. 6bd, respectively. The bursting trajectory (cyan) evolves in clockwise, as shown by the black arrow. Therefore, the burst is identified to terminate at the SH1 point, which exhibits three characteristics: (1) The burst terminates at the SH1 point (the blue star nearly overlaps with blue diamond); (2) The amplitude of the limit cycle (pink) is consistent with that of the burst; (3) The burst begins not at the SN1 point but a phase after the SN1 point (the red star is far away from the red diamond). Due to that the C lum is the slowest variable in the burst duration, the termination phase of the burst can be explained with the dynamics of the fast subsystem with C lum taken as the slow variable. This is the cause that the burst terminates at the SH1 point and the amplitude of the burst is consistent with that of the stable limit cycle. In addition, due to that C lum is not the slowest variable during the quiescent state, the quiescent state terminates not at the SN1 point (red diamond), i.e., the burst begins not from the SN1 point. Therefore, the dynamics of the termination phase of the quiescent state and initiation phase of the burst should be studied in the fast subsystem (V , n, C lum ) with C in taken as the bifurcation parameter in the following paragraphs, due to that C in is the slowest variable during the quiescent state.
3.5 Possible saddle-node bifurcation for the beginning phase of burst: manner 2 In the present subsection, C in is chosen as slow variable, and the equations of other three variables V , n, and C lum (Eqs. (1), (2) and (4)) combine to form the fast subsystem.
The difficulty introduced in Sect. 3.3 is encountered when the fast-slow variable dissection analysis is performed for the manner 2. Either of Eqs. (1) and (2) does not contain the variable C lum . Therefore, Eqs. (1) and (2) of V and n are independent of C lum , and Eq. (4) of C lum is independent of V or n. However, Eqs. (1) and (2) are dependent of C in , and Eq. (4) is dependent of C in . Therefore, Eqs. (1) and (2) are uncoupled to Eq. (4), although their dynamics are determined by the same parameter C in .
For the equilibrium point corresponding to a fixed C in of the fast subsystem, the equilibrium point of Eqs. (1) and (2) is acquired as (V (C in ), n(C in )) and of Eq. (4) is C lum (C in ), and then the equilibrium point of the fast subsystem for Eqs. (1), (2), and (4) is (V (C in ), n(C in ), C lum (C in )). With respect to C in , (V (C in ), n(C in ), C lum (C in )) forms the equilibrium curve (navy) of the fast subsystem (Eqs. (1), (2), and (4)), called equilibrium curve E2 (navy) in the present paper, as illustrated in Fig. 7. The projections of the curve E2 in the (C lum , C in , V ) and (C lum , C in , n) spaces are depicted in Fig. 7a, b, respectively, exhibiting Zshaped curves. A saddle-node bifurcation point appears as the middle branch (saddle) intersects with the lower branch (stable node), which is labeled as SN2 (red dot).
However, Eqs. (1) and (2) (V and n) exhibit a stable limit cycle, which is independent of Eq. (4) (C lum ), due to that Eqs. (1) and (2) and Eq. (4) are uncoupled. The projection of E2 (navy), maximal (LCU2, upper blue curve) and minimal (LCL2, lower blue curve) values of the limit cycle are depicted in Fig. 8a. A saddlehomoclinic bifurcation point emerges when the limit cycle intersects with the saddle (middle branch), which is labeled as SH2 (blue dot). The enlargement of Fig. 8a plotted with the trajectory of the period-11 bursting (cyan) forms Fig. 8b, and the projection on the (C in , n) plane is illustrated in Fig. 8c. In addition, the projection on the (C lum , C in ) space is depicted in Fig. 8d. Three characteristics can be found from Fig. 8b-d, which are introduced as follows: (1) The burst begins from the SN2 point in V , n, and C in directions (the red star and red dot are very close in Fig. 8b, c) instead of C lum direction (the red star is far from the red dot in C lum direction), which shows that the burst possibly begins and the quiescent state possibly terminates at the saddle-node bifurcation. The mismatch between the trajectory (navy) of the bursting and the SN2 point (red dot) of the fast subsystem in C lum direction is induced by the uncoupling relationship between Eqs. (1), (2) and Eq. (4). (2) The amplitude of the limit cycle, especially the minimal value (lower blue), is not consistent with (higher than) that of the burst, as illustrated in Fig. 8b, c. (3) The burst is identified to terminate not at the SH2 point (the blue star is very far from the blue dot).
In summary, the amplitude of the stable limit cycle is inconsistent with that of the burst, and the burst terminates not at the SH2 point, due to that C in is not the slowest variable in burst duration. The burst is identified to begin from the saddle-node bifurcation point in the V and n directions except for the C lum direction, due to that C in is the slowest variable in the quiescent state and the uncoupling characteristic between the differential equation of C lum and the differential equations of V and n. To acquire more comprehensive and deep  Fig. 9a to form Fig. 9b. Figure 9c represents the enlargement of Fig. 9b and the trajectory (cyan) of the period-11 bursting. Figure 9d is the enlargement of a part of Fig. 9c. As depicted in Fig. 9, the dynamics of the fast subsystem (V , n, C in ) are different from those of the fast subsystem (V , n, C lum ). For example, the curve E1 (red curve) is different from E2 (navy curve), SN1 (red diamond) is different from SN2 (red dot), and SH1 (blue diamond) is different from SH2 (blue dot). In c. E1 (red curve) and E2 (navy curve) represent the equilibrium point curves, SN1 (red diamond) and SN2 (red dot) represent the saddle-node bifurcations, and SH1 (blue diamond) and SH2 (blue dot) represent the saddle-homoclinic bifurcation points. LCU1 (upper pink curve) and LCU2 (upper blue curve) represent the maximal values of the limit cycles, and LCL1 (lower pink curve) and LCL2 (lower blue curve) represent the minimal values of the limit cycles. Blue and red stars represent the termination and initiation phases of the burst, respectively. (Color figure online) addition, the amplitudes (pink and blue) of the limit cycles are different. Furthermore, the SH1 (blue diamond) instead of SH2 is associated with the termination phase (blue star) of burst, and the amplitude (upper and lower pink) of the limit cycle for manner 1 instead of manner 2 (upper and lower blue) matches that of the burst (cyan). However, SN2 (red dot) is related to the beginning phase (red star) in the V , n, and C in directions except for the C lum direction, and SN1 (red diamond) is very far from the beginning phase of the burst (red star). Therefore, the bifurcation associated with the initiation phase of the burst is acquired by manner 3 of fast-low variable dissection method (C in and C lum taken as the slow variables) and the relationships of the dynamics between different fast subsystems are identified in the following paragraphs. 3.7 Saddle-node bifurcation for the beginning phase of the burst: manner 3 In this subsection, both C in and C lum are taken as the slow variables, and the fast subsystem is composed of Eqs. (1) and (2). The fast subsystem contains C in , but does not contain C lum . Therefore, the bifurcations of the fast subsystem are dependent of C in , but independent of C lum . For the fast subsystem, the bifurcations are depicted in Fig. 10. A Z-shaped curve of equilibrium (V , n) appears, called E in the present paper, as shown by the olive curve in Fig. 10a. A Hopf bifurcation point labeled as H (black circle) emerges on the upper branch. A stable limit cycle with maximal (LCU, green curve) and minimal (LCL, green curve) values emerges via the Hopf bifurcation point. When the limit cycle contacts with the middle branch, a saddle-homoclinic bifurca- tion called SH point (blue circle) emerges. The lower branch and the middle branch of the equilibrium curve intersect to form a saddle-node bifurcation labeled as SN (red circle). The projections of the curve E on the (C in , V ) plane and on the (C in , n) plane are depicted in Fig. 10b, c, respectively. The dynamics of the fast subsystem Eqs. (1) and (2) is independent of C lum , as shown in the (V , n, C in ) space in Fig. 10a. Therefore, the dynamics for the fast subsystem cannot be build a relationship to the trajectory (cyan) of period-11 bursting in the C lum direction. To overcome such a problem, the dynamics shown in Fig. 10 can be extended by the introduction of the dimension C lum . With changing C lum , (V , n, C in ) values along C lum -axe remain unchanged. For example, for the projection of the equilibrium point curve E (olive curve in Fig. 10b) on the (C in , V ) plane, after the introduction of dimension C lum and with (C in , V ) value along the C lum -axe unchanged, a surface in the (C lum , C in , V ) space forms, as shown in Fig. 11a. Such a surface extended from the curve E is called S here. Corresponding to the curve E which has the upper, middle, and lower branches, the surface S is also composed of three branches, called S U , S M , and S L , respectively. The surface is vertical to the (C in , V ) plane. Such a surface exhibits a relationship to the trajectory of the period-11 bursting in the (C lum , C in , V ) space.
Similarly, the saddle-node bifurcation point depicted in Fig. 10b can be extended to from a curve vertical to the plane (C in , V ), which is called SNc, as shown by the dashed yellow curve in Fig. 11a. Due to that the equations of the fast subsystem (V , n) do not contain C lum , the C lum value for the saddle-node bifurcation cannot be acquired in the fast subsystem (V , n). There- fore, the SNc curve exhibits fixed value of V , n, and C in along the C lum -axe.
Furthermore, the curves of the maximal (green curve) and minimal (green curve) values of the stable limit cycle in Fig. 10b can be extended to the surfaces upper and lower to the upper branch of the S surface S U , which is called L U and L L , respectively, as shown by the upper and lower olive surfaces in Fig. 11b. The intersection between the surface L L and the middle branch of S surface S M forms an SH curve (solid yellow curve), which is called SHc. The SHc curve (solid yellow curve) is the extension of the SH point along the C lum direction in Fig. 10. The enlargement of Fig. 11b around the SNc (dashed yellow curve) and SHc (solid yellow curve) points is shown in Fig. 11c. Furthermore, Fig. 11c around the middle and lower branches of S surface S M and S L is enlarged, as illustrated in Fig. 11d. These figures are helpful for understanding the dynamics of the period-11 bursting in the (C lum , C in , V ) space.
The relationship between the trajectory of the period-11 bursting and dynamics of the fast subsystem (V , n) in the (C lum , C in , V ) space can be acquired from Fig. 12. The trajectory (cyan) of the period-11 bursting and the dynamics such as the bifurcation curves and surfaces of the fast subsystem are illustrated in Fig. 12a. A part of Fig. 12a is enlarged and is depicted in Fig. 12b, and Fig. 12c is the back of Fig. 12b. After running across the middle branch S L at a certain point near the saddle-node curve SNc (dashed yellow curve), the quiescent state of the bursting terminates and the burst begins, as shown by the red star. Therefore, the period-11 bursting begins from the saddle-node bifurcation in the (C lum , C in ) space indeed, as shown in Fig. 12b. The result shows that the termination phase of the quiescent state, i.e., the initiation phase of the burst, is related to the saddle-node bifurcation. However, C lum values for the bifurcations cannot be acquired in the fast subsystem (V , n), due to that the equations of the fast sub- Fig. 12 The fast-slow variable dissection to the period-11 bursting (cyan) with C lum and C in taken as slow variables (manner 3). The surfaces of equilibrium and stable limit cycle and the bifurcation curves are the same as those in Fig. 11. a The framework of dynamics in the (C lum , C in , V ) space; b the enlargement of  Fig. b). The dashed yellow curve and the solid yellow curve represent the saddle-node curve SNc and the saddle-homoclinic curve SHc, respectively. The red and blue stars represent the initiation and termination phases of the burst, respectively. L U (upper green surface) and L L (lower green surface) represent the extended surfaces of the maximal and minimal values of the stable limit cycle, respectively. (Color figure online) system (V , n) do not contain the variable C lum , which is the cause that the burst begins from the SNc curve with a fixed value of V , n, and C in (exclusive the C lum value). In addition, the termination phase (blue star) of the burst is far from the SHc curve (solid yellow), which shows that the burst terminates not at the SHc curve. The amplitude of the limit cycle (L U and L L ) does not match that of the burst.
In summary, the burst of the period-11 bursting begins from the saddle-node bifurcation curve SNc of the fast subsystem (V , n). However, the burst terminates not at the SHc of the fast subsystem (V , n) or the SH2 of the fast subsystem (V , n, C lum ), but at the SH1 of the fast subsystem (V , n, C in ). In the following paragraphs, the comparisons of dynamics between different fast subsystems are presented.

Comparison between the dynamics of fast subsystems (V , n) and (V , n, C in )
The surface of equilibrium point S, the saddle-node bifurcation curve SNc (dashed yellow curve), and the saddle-homoclinic bifurcation curve SHc (solid yellow curve) for the fast subsystem (V , n) are shown in Fig. 13a. The equilibrium curve E1 and the saddlenode bifurcation point SN1 (red diamond) of the fast subsystem (V , n, C in ) are plotted with Fig. 13a to form Fig. 13b. The equilibrium point of the fast subsystem (V , n, C in ) is a subset of the corresponding bifurcation curve of the fast subsystem (V , n), and the saddle-node bifurcation point of the fast subsystem (V , n, C in ) is a subset of the saddle-node bifurcation curve SNc of the fast subsystem (V , n). Therefore, the equilibrium curve E1 is on the equilibrium surface S, and the SN1 point is on the SNc curve (dashed yellow curve). The cause for the relationships between the equilibrium points is given as follows. For the fast subsystem (V , n), the equilibrium points (V 0 , n 0 ) on the curve E (Fig. 10a) are the solution of the equations dV /dt = 0 and dn/dt = 0 at fixed C in values, due to that the fast subsystem (V , n) does not contain C lum , then, The points (V 0 , n 0 ) change with respect to C in to form the equilibrium point curve E shown in Fig. 10a, and V 0 values change with respect to C in to form the equilibrium curve depicted in Fig. 10b. The V 0 curve is extended along C lum direction to form the equilibrium surface S shown in Figs. 11a and 13a, and the SNc curve is formed by the extension of SN point in Fig. 10b along C lum direction. For the fast subsystem (V , n, C in ) with C lum regarded as the slow variable, the equilibrium points at different C lum values, i.e., the red curve E1 in Fig. 13b, are the solution of the equations dV /dt = 0, dn/dt = 0, and dC in /dt = 0, which are located on the equilibrium surface S, due to that the points on the red curve E1 are the solution of the equations dV /dt = 0 and dn/dt = 0 restricted by the equation dC in /dt = 0. The equilibrium curve (red) E1 of the fast subsystem (V , n, C in ) is related to C lum . Then, the equilibrium curve (red) E1 intersects with the SNc curve (independent of C lum ) to form the SN1 point.
However, unlike the equilibrium curve/surface, the limit cycle of the fast subsystem (V , n) is not directly related to the limit cycle of the fast system (V , n, C in ), due to that the two fast subsystems are different. In Fig. 13c, the extended surfaces of maximal (L U , upper green surface) and minimal (L L , lower green surface) Fig. 14 The comparison of fast-slow variable dissection between manner 1 (C lum taken as the slow variable) and manner 3 (both C in and C lum taken as the slow variables). a The trajectory (cyan) of the period-11 bursting, dynamics of the fast subsystem (V , n) including the equilibrium surface S (S U , S M , and S L ), the saddle-node bifurcation curve SNc (dashed yellow curve), and the saddle-homoclinic bifurcation curve SHc (solid yellow curve), and the dynamics of the fast subsystem (V , n, C in ) con-taining the limit cycle (LCU1 and LCL1, pink), the equilibrium curve E1 (red curve), the saddle-node bifurcation SN1 (red diamond), and the saddle-homoclinic bifurcation point SH1 (blue diamond) of the fast subsystem (V , n, C in ); b the enlargement of Fig. a; c the enlargement of Fig. b in a different angle; d Fig. c shown in a different angle (the back of Fig. c). The red and blue stars represent the initiation and termination phases of the burst, respectively. (Color figure online) values of the limit cycle of the fast subsystem (V , n) and the stable limit cycle with maximal (upper pink curve) and minimal (lower pink curve) values and the saddlehomoclinic bifurcation point SH1 (blue diamond) are superposed in Fig. 13b. Figure 13d is the enlargement of Fig. 13c, showing that the SH1 point (blue diamond) of the fast subsystem (V , n, C in ) is not located on the curve SHc (solid yellow curve) of the fast subsystem (V , n). In addition, the amplitude of the limit cycle of the fast subsystem (V , n, C in ), i.e., LCU1 and LCL1, is not consistent with that of the extended surfaces of the stable limit cycle L U and L L , due to that the L L surface of the fast subsystem (V , n) contains the SHc curve rather than LCL1 (minimal value of limit cycle) of the fast subsystem (V , n, C in ).
The trajectory (cyan) of the period-11 bursting is plotted with the dynamics of the fast subsystems (V , n) and (V , n, C in ), as depicted in Fig. 14a. Figure 14b represents the enlargement of Fig. 14a, Fig. 14 c is the enlargement around the SNc (dashed yellow curve) and SHc (solid yellow curve) of Fig. 14b, and Fig. 14d is the back of Fig. 14c. It can be found that the burst begins from the saddle-node bifurcation curve SNc (dashed yellow curve) of the fast subsystem (V , n), as illustrated in Fig. 14c. However, the initial phase of the burst (red star) and the SN1 point (red diamond) of the fast subsystem (V , n, C in ) have nearly equaling values of C in , which is not presented in Fig. 5, although they have different values of C lum . The burst begins not from the SN1 point (red diamond) of the fast subsystem (V , n, C lum ), due to that C in is not the slowest vari- able during the burst. As can be found from Fig. 14b, the maximal (upper pink curve) and minimal (lower pink curve) values of the burst (cyan) correspond to the maximal (upper pink curve) and minimal (upper pink curve) values of the limit cycle of the fast subsystem (V , n, C lum ), respectively. More importantly, the burst terminates at the SH1 point (blue diamond) of the fast subsystem (V , n, C lum ) instead of the SH curve SHc (solid yellow curve) of the fast subsystem (V , n), as depicted in Fig. 14b-d. 3.9 Comparison between the dynamics of the fast subsystem (V , n) and (V , n, C lum ) The equilibrium surface S, the saddle-node bifurcation curve SNc (dashed yellow curve), and the saddlehomoclinic bifurcation curve SHc (bold yellow curve) of the fast subsystem (V , n) are shown in Fig. 15a. Figure 15b represents Fig. 15a plotted with the equilibrium curve E2 (navy) and saddle-node bifurcation point SN2 (red dot) of the fast subsystem (V , n, C lum ). Similarly to Fig. 13, the curve E2 (navy) is on the surface S and the SN2 point (red dot) is on the SNc curve (dashed yellow curve), and the cause for the relationships of the equilibrium point/curve/surface between the two fast subsystems resembles that of Fig. 13. However, unlike the equilibrium curve/surface, the limit cycle of the fast subsystem (V , n) has no direct relations to that of the fast subsystem (V , n, C lum ). In Fig. 15c, the stable limit cycle with maximal (upper blue curve) and minimal (lower blue curve) values and the saddle-homoclinic bifurcation point SH2 (blue dot) of the fast subsystem (V , n, C lum ), and the extended surfaces of maximal value (L U , upper green surface) Fig. 16 The comparison of fast-slow variable dissection between manner 2 (C in taken as the slow variable) and manner 3 (both C in and C lum taken as the slow variables). a The trajectory (cyan) of the period-11 bursting, the dynamics of fast subsystem (V , n) including the equilibrium surface S (S U , S M , and S L ), the saddle-node bifurcation curve SNc (dashed yellow curve), and the saddle-homoclinic bifurcation curve SHc (solid yellow curve), and the dynamics of the fast subsystem (V , n, C lum ) con-taining the equilibrium curve E2 (navy curve), stable limit cycle (LCU2 and LCL2, blue), saddle-node bifurcation point SN2 (red dot), saddle-homoclinic bifurcation point SH2 (blue dot) of the fast subsystem (V , n, C lum ); b the enlargement of Fig. a;  and minimal value (L L , lower green surface) of the stable limit cycle of the fast subsystem (V , n), are added to Fig. 15b. Figure 15d represents the enlargement of a part of Fig. 15c. The SH2 point (blue dot) of the fast subsystem (V , n, C lum ) is on the curve SHc (solid yellow) of the fast subsystem (V , n), due to that the curve is extended from the SH2 point. And the LCU2 and LCL2 curves are on the L U and L L branches of the surface S, respectively, due to that the surfaces are extended from the corresponding curves.
The trajectory (cyan) of the period-11 bursting is plotted with the dynamics of the fast subsystems (V , n) and (V , n, C lum ), as depicted in Fig. 16. Figure 16b is the enlargement of Fig. 16a, Fig. 16c represents the further enlargement of Fig. 16b around the SNc (dashed yellow curve) and SHc (solid yellow curve) curves, and Fig. 16d is the back of Fig. 16c. The burst begins from the saddle-node bifurcation curve SNc (dashed yellow curve) of the fast subsystem (V , n), as shown in Fig. 16c. The initial point of the burst (red star) and the SN2 point (red dot) of the fast subsystem (V , n, C lum ) have nearly equaling values of C in but different values of C lum . The burst terminates not at the SH2 point (blue dot) of the fast subsystem (V , n, C lum ) or the saddlehomoclinic bifurcation curve SHc (solid yellow curve) of the fast subsystem (V , n), as illustrated in Fig. 16b, d. As depicted in Fig. 16b, the maximal and minimal values of the burst (cyan) do not correspond to those of the limit cycle (upper and lower blue curves) of the fast subsystem (V , n, C lum ). All these results further show that the termination of the burst cannot be explained with the fast subsystems (V , n) or (V , n, C lum ).

Fig. 17
The comparisons between 3 manners of fast-slow variable dissection to the period-11 bursting. a The trajectory of the period-11 bursting (cyan), the beginning phase (red star) and ending phase (blue star) of the burst, and the dynamics of the fast subsystem (V , n, C in ) (equilibrium point curve E1 (red curve), saddle-node bifurcation point SN1 (red diamond), saddle-homoclinic bifurcation point SH1 (blue diamond), minimal value of the limit cycle LCL1 (lower pink curve)), the dynamics of the fast subsystem (V , n, C lum ) (equilibrium point curve E2 (navy curve), saddle-node bifurcation point SN2 (red dot), saddle-homoclinic bifurcation point SH2 (blue dot), minimal value of the limit cycle LCL2 (lower blue curve)), and the dynamics of the fast subsystem (V , n) (middle and lower branches of the extended equilibrium point surface S M (olive surface) and S L (olive surface), saddle-node bifurcation curve SNc (dashed yellow curve), saddle-homoclinic bifurcation curve SHc (solid yellow curve), and minimal value of the extended surface of the limit cycle L L (lower green surface)); b Fig. a shown in a different angle; c The bifurcations related to the beginning (red star) and ending (blue star) phases of the burst in Fig. a;  The surfaces, bifurcation curves, bifurcation points of three fast subsystems, and the trajectory of the period-11 bursting around the saddle-node bifurcation curve SNc are depicted in Fig. 17a. Figure 17b represents the back of Fig. 17a. The largest amplitudes of the burst are not shown here due to that they are introduced clearly in the previous subsections. The relationships between these surfaces, curves, and points mentioned above can be clearly found from Fig. 17, and the dynamics of 3 fast subsystems are also listed in Table 1. Figure 17c, d represents the dynamics related to the begging and ending phases of the burst from front and back sides, respectively. There are two most important points about the dynamics of the period-11 bursting, which are described as follows: (1) The burst begins from a certain point on the saddlenode curve SNc (dashed yellow curve) of the fast subsystem (V , n), as illustrated in Fig. 17a. The C lum of the point cannot be acquired from the fast subsystem (V , n) due to that the fast subsystem (V , n) does not contain C lum , which is an important characteristic of the model containing C in and C lum . The C lum influences the dynamics of the bursting via the modulations of C lum to C in . The (C in , V , n) value of the termination point of the burst is determined by the fast subsystem. Correspondingly, the burst begins not at the saddle-node bifurcation SN2 of the fast subsystem (V , n, C lum ) due to the uncoupling characteristic between the equation of C lum and equations of V and n. In addition, the burst starts not at the saddle-node bifurcation SN1 of the fast subsystem (V , n, C in ), due to that C lum is not the slowest variable during the quiescent state. The results show that the initial phase of the burst can be explained with the fast subsystem (V , n) through the fast-slow variable dissection analysis with two slow variables C lum and C in related to the calcium concentrations.
(2) The burst terminates at the saddle-homoclinic bifurcation point SH1 (blue diamond) of the fast subsystem (V , n, C in ), and the amplitude of the limit cycle (pink) matches that of the burst, as depicted in Fig. 17c, d, which shows that the burst dynamics can be explained with the dynamics of the fast subsystem (V , n, C in ), due to that C lum is the slowest variable during the burst. Therefore, when C in is taken as the slow variable or both C in and C lum are chosen as the slow variables, the burst dynamics cannot be well explained, which is the cause that the burst terminates not at the saddlehomoclinic bifurcation point SH2 of the fast subsystem (V , n, C lum ) or the saddle-homoclinic bifurcation curve SHc of the fast subsystem (V , n).
In the present paper, the period-11 bursting is chosen as the representative to show the results. Such results for bursting are independent of parameter values of v c . In fact, if period-4 bursting for v c = 110 mV (Fig. 1a) is chosen, the results are similar. To avoid possible repetitions, the results for the period-4 bursting are not addressed in the present paper.

Conclusion
Bursting is an important nonlinear behavior of neurons related to physiological or pathological functions. Therefore, identification of nonlinear dynamics underlying bursting is a very important issue [1][2][3][4][5][6][60][61][62]. The dynamics of bursting modulated by one slow variable has been well explained with the bifurcations acquired by the fast-slow variable dissection in the previous studies [1,22,[27][28][29]. At present, identification of bursting dynamics modulated by two slow variables is still an open problem. Two cases of bursting modulated by two slow variables have been studied with fast-slow variable dissection method in the recent studies [54,58]. In the present paper, the dynamics of the third case of bursting modulated by two slow variables related to the calcium concentrations, C in and C lum , is acquired in a four-dimensional model to describe the dynamics of a neural pacemaker related to the pathological pain. The results exhibit significance in the following 4 aspects.
Firstly, a feasible process to perform fast-slow variable dissection is proposed to analyze the dynamics of bursting modulated by two slow variables C in and C lum . When traditional process of fast-slow variable dissection method with two slow variables is performed to analyze the case-3 bursting, the problems are encountered due to that the fast subsystem (V , n) does not contain C lum . The problems are overcome by introducing the dimension C lum . For the fast subsystem (V , n), we can only acquire one-parameter bifurcations with respect to C in , which is independent of the parameter C lum . Based on the independence, the one-parameter bifurcation points can be extended along the dimension C lum to form bifurcation curves in two-parameter space (C lum , C in ). Although such extended bifurcations in two-parameter space are independent to C lum , which can be easily used to identify the bifurcations corresponding to the initiation phase of the burst for case-3 bursting.
Secondly, such a method can be extensively employed to analyze bursting patterns modulated by C in and C lum . The results for the bursting (case-3) in the present paper are acquired in the four-dimensional model of pacemaker with C in and C lum . In fact, for other nervous system or neuronal model such as pre-Bötzinger (PB) complex and β-cells [23,36,53,54], there are two slow variables C in and C lum . The method used in the present paper can be extensively used to analyze the bursting of these models. For these models, except for the existence of two slow variables C in and C lum , the equation to describe dV dt contains C in instead of C lum . However, such a method is different from the fast-slow variable dissection method with two slow variables used in Ref. [54] (case-1) and Ref. [58] (case-2), due to that the characteristics of the slow variables or the relationship between the slow variables are different.
Thirdly, the initial phase of the burst for case-3 bursting is identified with the novel method of fast-slow variable dissection. The burst is identified to begin from the saddle-node bifurcation curve of the twodimensional fast subsystem (V , n). However, the termination phase mismatches the saddle-homoclinic bifurcation curve of the subsystem (V , n), due to that C in is much slower than C lum in the quiescent state.
Last, the burst is identified to terminate at the saddlehomoclinic bifurcation point of the fast subsystem (V , n, C in ) with C lum taken as the bifurcation parameter, due to that C lum is the slowest variable in the burst duration. On the one hand, the present paper shows that the combination of 3 manners of fast-slow variable dissection is more effective than one out of the three manners. On the other hand, the deep understandings for the complex bursting are helpful for the modulation to the bursting behaviors related to the pathological pain.
Data availability The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.