Timescales and Mechanisms of Sigh-Like Bursting and Spiking in Models of Rhythmic Respiratory Neurons

Neural networks generate a variety of rhythmic activity patterns, often involving different timescales. One example arises in the respiratory network in the pre-Bötzinger complex of the mammalian brainstem, which can generate the eupneic rhythm associated with normal respiration as well as recurrent low-frequency, large-amplitude bursts associated with sighing. Two competing hypotheses have been proposed to explain sigh generation: the recruitment of a neuronal population distinct from the eupneic rhythm-generating subpopulation or the reconfiguration of activity within a single population. Here, we consider two recent computational models, one of which represents each of the hypotheses. We use methods of dynamical systems theory, such as fast-slow decomposition, averaging, and bifurcation analysis, to understand the multiple-timescale mechanisms underlying sigh generation in each model. In the course of our analysis, we discover that a third timescale is required to generate sighs in both models. Furthermore, we identify the similarities of the underlying mechanisms in the two models and the aspects in which they differ.


Introduction
Many years of experimental work have elucidated a range of properties of the neuronal circuits involved in breathing. While there is clear evidence that one or more neuronal populations in the mammalian brain stem, including the well-studied pre-Bötzinger complex (pre-BötC) [1], are capable of generating respiratory-related rhythms, there is a lack of consensus about the transfer from rhythm generation to patterned motor output. While some have presented evidence that modulations of the interactions within a common rhythmogenic core can account for a variety of patterns of activity observed in recordings from respiratory pathways [2], others have argued that rhythm generation and pattern generation are distinct, with one or more pattern-generating respiratory populations shaping stereotyped baseline rhythms into different outputs [3,4]. Recent computational work has shown, for example, that a distinct pattern generator is not necessary to explain the observation of mixed-mode oscillations composed of intermingled bursts and burstlets in respiratory neuron population readout nerves [5]. Nonetheless, further study is needed to help distinguish between these competing views.
One familiar example of a non-standard respiratory output pattern is the sigh. Sighs can be distinguished in neuronal activity; indeed, recordings of population activity from the ventral respiratory group (VRG) identified two distinct patterns of inspiratory activity under conditions of normal oxygenation: fictive eupnea (typical respiration) and low-frequency fictive sighs. As described by [6], fictive sighs occur periodically in the in vitro transverse medullary slice preparation from mice containing the pre-BötC within the VRG. Each fictive sigh consists of a biphasic activity burst that is larger in amplitude, longer in duration, and occurs at a lower frequency than eupneic bursts. Despite the experimental accessibility of sighs, the rhythmogenic mechanisms underlying sigh generation remain largely unknown [6,7]. The biphasic aspect of the sigh, with an initial phase that is identical to a normal eupneic burst and a later high-amplitude phase, could result from the recruitment of a neuronal population distinct from the eupneic rhythm-generating circuit, or it could simply emerge due to the complex interplay of multiple-timescale processes within the core rhythm-generating circuit itself. The main goal of this paper is to use mathematical tools to elucidate the multiple-timescale mechanisms underlying sigh generation in two recent computational models [8,9], one of which represents each of the competing hypotheses about pattern generation. In doing so, we will highlight the ways in which the dynamic mechanisms in the two models are in fact similar as well as ways that they can be distinguished in future experimental studies, to help determine whether or not separate pattern generating components complement rhythm generators in producing respiratory outputs.
Example sigh patterns produced by the two models appear in Fig. 1. The model yielding the solution shown in Fig. 1A includes fast-spiking currents in addition to rhythmic burst generation, while the one associated with Fig. 1B does not, leading to the significant quantitative differences between their outputs. For convenience, we refer to such solutions as sigh-like bursting (SB, Fig. 1A) and sigh-like spiking (SS, Fig. 1B), respectively. Of course, the meaning of a 'spike' is quite different across the two models, representing a single action potential in the SB case and an entire active  [8] and [9], respectively. Patterns repeat periodically. In (B), we only plot a few of the many eupneic cycles occurring between sighs, to obtain better resolution in displaying the sigh period in the SS model. Nonetheless, both feature high-amplitude, low-frequency, long-duration events emerging periodically on the top of higher-frequency baseline patterns. Indeed, a comparison of the patterns in Fig. 1 suggests that the fundamental dynamic mechanisms underlying both may be similar, analogously to the comparison between square-wave (or fold-homoclinic) bursting with spikes and relaxation oscillations without them. One of the contributions of our analysis will be to determine the extent to which this analogy holds, which will clarify the relation of these models for subsequent studies.
Experimental studies in rodent medullary slices containing the pre-BötC have identified two biophysical mechanisms that could potentially contribute to the generation of rhythmic bursting, one based on the persistent Na + current (I NaP ), and the other involving the voltage-gated Ca 2+ current (I Ca ) and the Ca 2+ -activated nonspecific cation current (I CAN ), activated by intracellular Ca 2+ , which may be accumulated from a variety of sources [10,11]. Past computational work showed that the interactions of these burst mechanisms could yield a form of mixed bursting (MB) output with significant qualitative similarity to the SB pattern shown in Fig. 1A [12][13][14]. By applying methods of dynamical systems theory to a singlecompartment model of a pre-BötC inspiratory neuron in [15], we explained in full detail how the MB solution results from these currents. In this paper, we generalize the analysis and methods for studying MB solutions in [15] to the more complicated model SB and SS models presented in Sect. 2 in order to uncover the mechanisms underlying the these patterns. Our approach is geometric and is based on studying reduced subsystems of an original model that evolve on particular timescales. It is not rigorous, in that we assume that there are abrupt transitions between timescales and we do not prove any results; moreover, we will sometimes make approximations, such as treating a nullcline that only weakly depends on a parameter as fixed under variations of that parameter. Nonetheless, this approach has a long history of providing powerful insights, for example, in work ensuing from the classical dissection of minimal bursting models by [16] and in the study of coupled neuronal oscillators and bursters (e.g, [17][18][19][20] and many others since). An important aspect of this dynamical systems approach is that it can decompose a solution pattern into a sequence of dynamic epochs evolving on different timescales and bifurcations of subsystems that underlie transitions between these epochs [21]. Thus, the approach distills out a set of key features that can be used to classify solutions and to objectively compare solutions that differ quantitatively and come from models that are superficially quite different, as in the famous classification of neuronal bursters [16,22].
While the MB solution appears to involve at least three timescales based on its time course, we obtained a non-intuitive result in [15] that the core mechanism underlying the robust production of the MB pattern is in fact an interaction of only two timescales. In this paper, after presenting the SB and SS models in Sect. 2, we similarly analyze the timescale-interaction mechanisms that support SB solutions (in Sect. 3) and SS solutions (in Sect. 4). In both cases, we determine that timescales must remain separated into three distinct classes for the solutions of interest to arise. In Sect. 5, we compare these two sighing solutions to highlight the similarities and differences in the mechanisms and timescale interactions involved in producing them, and we conclude with a discussion in Sect. 6.

Sigh-Like Bursting Model
We consider two recent models for neurons in the pre-BötC. Jasinski et al. [8] presented a relatively detailed model in the Hodgkin-Huxley framework for pre-BötC neurons and showed that a synaptically coupled population of these neurons, with heterogeneous parameter values, can generate SB solutions (Fig. 1A), whereas a single model pre-BötC neuron without synaptic inputs cannot produce a sighing rhythm. Beyond issues of synaptic coupling, the model by Jasinski et al. [8] is more complicated than the MB model [15] in two important ways. First, in addition to I NaP and I CAN , the Na + /K + pump plays an important role in the generation of activity patterns in this model. Furthermore, we have considered only one-directional coupling, from calcium concentration to the voltage dynamics, in [15], such that the MB model can be thought of as one oscillator forcing another (see also [23]). In the model developed by Jasinski et al., however, the membrane potential and the cytoplasmic Ca 2+ concentration can each influence the evolution of the other.
To facilitate the identification of the essential mechanisms underlying the SB behavior and the assessment of how to group the timescales involved, we consider a synaptically self-coupled single-compartment model neuron based on the model presented in [8], which we refer to as the Jasinski model: The neuronal membrane potential (V ) is governed by a set of membrane ionic currents, as shown in (1a). C is neuronal membrane capacitance, which we set to 36 pF, and t is time. The ionic currents in the model include fast Na + current (I Na ), persistent Na + current (I NaP ), delayed rectifier K + current (I K ), high-voltage-activated Ca 2+ current (I Ca ), Ca 2+ -activated nonspecific cation current (I CAN ), Na + /K + pump current (I Pump ), leak current (I L ) and excitatory synaptic current (I SynE ), which combines self-coupling with tonic drive from respiratory feedback regions; see Table 1.
Activation (m) and inactivation (h) variables for most ionic channels are governed by Eq. (1b), where steady-state activation (m ∞ ) and inactivation (h ∞ ) functions and time constants are described as in Table 2. Notice from the bottom row in Table 2 that unlike the other currents, I CAN activates instantaneously. This activation (m CAN ) depends on the intracellular calcium concentration (Ca i ) and is voltage-independent. The parameters for these currents are specified in Table 3.
The dynamics of the total intracellular Ca 2+ concentration within the cell (Ca tot ) and intracellular concentration of free Ca 2+ (Ca i ) are described by (1c) and (1d), respectively, and these are intimately linked with l, the fraction of IP 3 receptors that are not inactivated by calcium, which is governed by (1e). Calcium dynamics is influenced by voltage through the first term in the right-hand sides of (1c) and (1d), −α Ca I Ca , which represents Ca 2+ influx from the extracellular space through voltagegated Ca 2+ channels. In (1d), J ERin represents the flux of Ca 2+ per unit volume from the endoplasmic reticulum (ER) into the cytoplasm, which depends on l, and J ERout Table 2 Functions associated with activation and inactivation variables for system (1a)-(1g). We use the variable y when an expression corresponds to a set of variables Gating variables

Steady-state activation and inactivation
Time constants The values of parameters associated with the equations for Ca tot , Ca i and l are given by α Ca = 2.5 × 10 −5 mM/fC, τ Ca = 500 ms, K Ca = 2.5 × 10 −5 , L IP 3 R = 0.37/ms, P IP 3 = 31,000/ms, [13]. Additional description of model components has been given previously [13,14]. For convenience, we henceforth omit units from the model parameters and variables.

Sigh-Like Spiking Model
Toporikova et al. [9] recently designed a computational model for inspiratory pre-BötC neurons, based on an earlier model [13] that includes two different bursting mechanisms depending on I NaP and intracellular Ca 2+ , respectively. In contrast to the Jasinski model, this model does not include fast-spiking components (i.e., I Na and I K ) such that instead of generating SB solutions, it produces SS patterns. Under different parameter choices emphasizing distinct burst-generating mechanisms, this model generates oscillations following mainly the kinetics of I NaP or an even lowerfrequency sigh-like rhythm resulting mainly from slow Ca 2+ oscillations. Following the terminology used in [9], we can instantiate two copies of the model, one with each type of parameter set, and refer to the model equipped with parameters that support fictive sigh activity as the sigh compartment and the model with parameters that support eupneic activity as the eupnea compartment.
With an inhibitory synapse from the eupnea compartment to the sigh compartment and an excitatory synapse from sigh to eupnea, the output of the coupled model is the SS pattern (see Fig. 1B) composed of a large spike emerging periodically on the top of regular spikes, analogous to the long bursts separated by short bursts in the SB solution, which is consistent with experimental data obtained in vitro [6,7,24]. In fact, the role of the coupling between compartments is simply to coordinate the timing of the eupneic and sigh-like spikes. For example, the oscillation pattern occurring after removal of the coupling from the eupnea compartment to the sigh compartment appears in Fig. 2A, while the spike of the sigh compartment alone, without input, is shown on a different timescale in Fig. 2B. Furthermore, the generation of eupneic spikes can be understood in the same way as the regular bursts in the SB solution, which we consider in Sect. 3.1.1. Hence we focus entirely on the sigh compartment on its own in our analysis in Sect. 4. For simplicity, we will still denote the solution shown in Fig. 2B as the SS solution and we refer to the sigh compartment model as the Toporikova model, described by the following equations: with In this system, the terms J ER IN and J ER OUT appearing in the Ca i Eq. (3d) are the same as those used in the Jasinski model as given by Eqs. (2a)-(2c). Parameters related to them have been specified in Sect. 2.1, except [IP 3 ] = 1 μM. In Eqs. (3c)-(3d), λ is the ratio of ER to plasma membrane surfaces, which we set to 0.1. C m , the neuronal membrane capacitance, is 21 pF. Other parameter values related to currents and corresponding units for the Toporikova model are given in Table 4.

Sigh-Like Bursting in a Self-Coupled Pre-BötC Neuron
At the single-neuron level, with the self-coupling removed by setting g SynE = 0, model (1a)-(1g) can produce different intrinsic bursting patterns, depending on chosen parameter values. One type of bursting is based on the slow voltage-dependent inactivation of I NaP (as represented by the h NaP variable in the equation for I NaP in Table 1), whereas another type relies on intracellular Ca 2+ and depends on I CAN . There are several distinct burst-terminating mechanisms that can contribute, based, respectively, on the slow inactivation of I NaP , the activity-dependent accumulation of Na + followed by the action of the [Na + ] in -activated I Pump , and the Ca 2+ -dependent inactivation of IP 3 receptors.
With synaptic coupling between two or more neurons (i.e., g SynE > 0), the coupled cells are able to generate SB solutions. We consider the special case of a single selfcoupled cell, given by model (1a)-(1g), as a reduction of a coupled network. Similarly to findings in [8], numerical simulations of (1a)-(1g) indicate that assigningḡ Ca = 0 orḡ CAN = 0 does not affect regular bursting, but fully removes the sigh-like oscillations and hence the SB pattern. These effects suggest that the generation of sigh-like bursts in this model is I Ca /I CAN -dependent, while regular bursting is independent of I CAN . Furthermore, setting I Pump = 0 eliminates both regular bursting and sigh-like oscillations, implying that the Na + /K + pump is also critically involved. On the other hand, simulations withḡ NaP = 0 show that the full SB solution survives intact without the need for I NaP . The left panel of Fig. 3 demonstrates the evolution of several variables, Ca i , Na i , Ca tot and l, during an SB cycle. Note that there exist low-amplitude Ca i and Na i transients during regular bursts, eventually followed by an abrupt increase of intracellular Ca 2+ , which Na i appears to follow more slowly. Meanwhile, Ca tot and l, which only interact directly with each other and Ca i , both accumulate until Ca i jumps up and then start decreasing. The right panel of Fig. 3 shows in a magnified view of a regular burst that Ca i , Na i and Ca tot engage in small oscillations during the spiking phase, consistent with the fact that these variables receive input from V ; they tend to increase throughout this phase and to decrease while V is not spiking.
We aim to understand the mechanisms underlying this SB solution and to elucidate the timescales that are needed to produce it. The bidirectional coupling between membrane potential and cytoplasmic Ca 2+ concentration, as well as the highdimensionality of model (1a)-(1g), are complications not present in previous related analyses [15,23]. To achieve these goals, we will nondimensionalize Eqs. (1a)-(1g) to reveal the presence of different timescales; determine how to group timescales that are present; implement geometric singular perturbation theory (GSPT) to set up reduced systems based on separation of timescales; and use the reduced systems to explain the mechanisms underlying the dynamics of the SB solutions. In contrast to the related MB model [15], it will turn out that use of the averaging method will play a role in the analysis. By uncovering the mechanisms underlying the SB solution in the Jasinski model, we will conclude that, unlike the case with the MB solution studied previously, a third timescale actually is required to generate SB solutions in this model.

Analysis of Sigh-Like Bursting
Since the SB solution dynamics in the Jasinski model persists without I NaP , we can reduce (1a)-(1g) to an 11-dimensional system by removing m NaP and h NaP and setting I NaP = 0 in Eq. (1a). Henceforth, we still refer to the new lower-dimensional model as the Jasinski model. The fact that regular bursting depends on I Pump but not I CAN and I Ca suggests that we can decouple the I Pump -based burster from the calcium dynamics during the regular bursting phase of the SB solution. Hence, we can think of the Jasinski model as consisting of two subsystems: (V , y, Na i , s) (denoted as the voltage compartment) as a potentially bursting subsystem based on I Pump , and (Ca i , Ca tot , l) (denoted as the calcium compartment) as a potentially oscillating subsystem. We first investigate the I Pump -based mechanism underlying the regular burst generated by the voltage compartment through a bifurcation analysis and then study the effect of the calcium compartment on the resulting bifurcation structures to understand how the transition from regular bursts to the long burst happens. Of course, we will have to take into account the coupling from the voltage compartment to the calcium compartment to complete the analysis.
Our methods for analyzing the model will depend heavily on exploiting the presence of different time scales. As a first step, it is helpful to rescale the variables so that the important timescales can be explicitly identified and used to group variables into timescale classes that are different from the voltage and calcium compartment groupings, which are based on coupling structure and biology. To this end, we define new dimensionless variables (v, ca i , c tot , na i , τ ), and voltage, calcium, sodium and timescales Q v , Q c , Q na , and Q t , respectively, such that Note that y, s and l are already dimensionless in (1a)-(1g). Details of the nondimensionalization procedure, including the determination of appropriate values for Q v , Q c , Q na and Q t , are given in Appendix 1. From this process, we obtain a dimensionless system of the form with coefficients of derivatives on the left-hand sides as well as functions on the right-hand sides specified in Eqs. (14a)-(14g), and timescales for all variables shown in Table 5, both of which appear in Appendix 1. While v, gating variables m Na , h Na , m Ca , h Ca , m K , and s do not operate on exactly the same timescale quantitatively, it is clear that they are all relatively faster than the other variables. Hence we choose to group all of them as fast variables, to consider na i and ca i as slow, and to classify l and c tot as evolving on a superslow timescale. For simplicity, we abuse notation to now let y ∈ R 6 denote all the fast gating variables along with s. For each group of variables we can define a corresponding subsystem of equations with slower variables kept as parameters, as we have done in [23] and many others have done previously. We can also define a fast-slow subsystem of fast and slow variables together, and we can define separate fast and slow subsystems for the voltage compartment, since it includes slow na i . The bifurcation diagram for the fast subsystem of the voltage compartment, comprising variables (v, y, s) and decoupled from ca i by settingḡ Ca =ḡ CAN = 0, with the slow variable na i treated as a bifurcation parameter, is shown in Fig. 4A. It includes an S-shaped curve of equilibria (S) and a family of stable periodic orbits (P ) that initiates in a supercritical Andronov-Hopf (AH) bifurcation and terminates in a homoclinic (HC) bifurcation involving the middle branch of S as na i is increased. Hence, in the absence of calcium dynamics, this subsystem is capable of generating a square-wave bursting solution, which terminates via the accumulation of na i and subsequent activation of the Na + /K + pump. As part of our analysis of SB dynamics, we will in Sect. 3.1.1 consider what happens to this bursting, corresponding to the small bursts in the SB solution, once coupling from the calcium compartment to the voltage compartment is restored.
In the calcium compartment, the dynamics of ca i depends on the neuronal membrane potential v. We can represent this v-dependence by considering a family of ca i nullsurfaces, each defined for v fixed. While ca i is low, the projection of the trajectory to (ca i , c tot , l)-space exhibits small oscillations during each regular burst (see Fig. 4C and the calcium trace in Fig. 3). These oscillations correspond to the projected trajectory trying to move back and forth between the left branches of two extreme nullsurfaces as v oscillates between its minimum and maximum during the spiking phase of each burst ( Fig. 4B and C); the trajectory cannot make it all the way to the v max surface because the dynamcs of ca i is slower than that of v. As for the right branches of these two extreme surfaces, Fig. 4B shows that they lie close together, which results because I Ca depends only weakly on v for calcium large. As a result, if ca i is elevated, then the projected trajectory is constrained quite tightly between the two right branches of these nullsurfaces. We can observe that at the end of a cycle of the SB solution, the sigh-like burst is completed as the trajectory passes the curve of lower folds of the family of calcium nullsurfaces and jumps back to the left. What remains unclear about this loop is what bifurcation induces the jump-up of calcium, the understanding of which is crucial in illustrating the transition from regular bursts to the high-amplitude sigh-like burst. We consider this issue in Sect. 3.1.2, after first completing some additional analysis of the regular bursting phase with coupling from the calcium compartment to the voltage compartment restored.

Mechanisms Underlying Regular Bursting
Settingḡ Ca = 0.00065 andḡ CAN = 3 as given in Table 3 restores the coupling from calcium to voltage and yields an SB solution. An example of the coupling effect on the voltage compartment can be seen in Fig. 5A: an increase of ca i shifts the na i -dependent fast subsystem equilibria S to the right. In Fig. 5B, we project the first regular burst solution and the bifurcation diagram of the fast subsystem for ca i fixed at 8e−3, corresponding to its value at the beginning of this first small burst, onto (na i , v)-space. Also shown is the green (resp. blue) dashed line representing the na i values at which the homoclinic (resp. lower fold of equilibria) bifurcation occurs. Starting from the yellow star, the trajectory moves on the slow timescale associated with na i along the stable lower branch of S until it reaches the lower fold. After that, the trajectory jumps up to the stable periodic orbit branch and then moves to the right, since the trajectory stays above the na i -nullcline. Sometime after it crosses the homoclinic bifurcation at the na i value indicated by the green dashed line, the trajectory will jump down to the lower branch of equilibria, completing a small burst. This is essentially a square-wave burst, but notice that several more spikes occur after the green dashed line is passed. These spikes arise because during The curve of saddle-node bifurcations corresponding to the lower fold of the bifurcation diagram (blue), homoclinic bifurcation curve (green) and part of the trajectory (black) generated by (5a)-(5g) in (na i , ca i )-space. The HC curve splits the (na i , ca i )-space into two regions labeled as 'Active' and 'Silent', respectively. The part of the trajectory corresponding to the first burst, as shown in (B), is magenta. (D) A zoomed-in and enlarged view of (C) the first regular burst period, ca i progressively increases on a slow timescale; as a result, the bifurcation diagram also moves rightward on a slow timescale associated with the increase of ca i (Fig. 5A). Hence at the end of the burst, the homoclinic bifurcation actually occurs at some larger na i value to the right of the green dashed line in Fig. 5B, yielding several more spikes after the green dashed line. Such squarewave bursting solutions will repeat roughly until ca i starts to jump up to larger values as indicated in Figs. 3, 4. Understanding the persistence of the regular bursts and the mechanism by which a transition to the sigh-like burst occurs requires us to consider the effect of ca i on the voltage compartment. To do so, we use ca i as the second bifurcation parameter and allow both ca i and na i to vary in order to find the two-parameter bifurcation curves of the fast subsystem (v, y) in the (na i , ca i ) parameter plane that unify the results in Fig. 5A, as illustrated in Fig. 5C. The blue (resp. green) curve in this plane is the curve of lower fold (LF) (resp. homoclinic (HC)) bifurcations, which initiates (resp. terminates) each burst, as noted previously. Since the increase of ca i moves the bifurcation diagram to the direction of increasing na i , both the LF and the HC curves are positively sloped in (na i , ca i )-space. Within the same projection, the trajectory evolves leftward from the yellow star and it starts oscillating as it passes the LF curve (see Fig. 5C). These oscillations terminate when the trajectory reaches the HC bifurcation curve, which completes the first regular burst. Similarly, a sequence of subsequent regular bursts occurs, with the local maximum of na i progressively increasing due to the rightward drift of LF as ca i accumulates. The fact that the trajectory in (na i , ca i )-space crosses the LF and HC 15 times corresponds to the existence of 15 regular bursts between sighs (see Fig. 3). After these, bursting solutions give way to continuous spiking.
Based on the fast voltage compartment bifurcation structures in this section, we have seen that regular bursts occur as the slow variables ca i and na i traverse the phase space back and forth between the LF and HC curves. The reason why regular bursts switch to the sigh-like burst, however, has not been addressed. To figure this out, we notice that after multiple crossings of the HC curve and returns to quiescence in Fig. 5C, the trajectory projected to (na i , ca i ) space starts oscillating near the HC curve, instead of going back again to the quiescent state (see Fig. 5D for an enlarged view of oscillations near the HC curve). Moreover, this transition happens before ca i jumps up. Hence, the switch from regular bursts to the long sigh-like burst in the full system seems to correspond to the transition from bursting to tonic spiking in the fast-slow subsystem (v, y, na i , ca i ), rather than the jumping up of ca i in the calcium compartment as in the MB model [15]. We next consider the mechanism responsible for this transition.

Mechanisms Underlying the Transition from Regular Bursts to the Sigh-Like Burst
In the analysis up to this point, the superslow variables c tot and l have not yet been considered. It is natural to expect that the superslow evolution of these two variables may contribute to the switch from bursting to tonic spiking. A numerical simulation of the fast-slow subsystem over a range of c tot and l values ( Fig. 6A-C) suggests that superslow variables do play an important role in inducing a transition from bursting to tonic spiking in the fast-slow subsystem. In Fig. 6A-C, representative time traces for For small c tot , the fast-slow system is in a bursting state (Fig. 6A). When c tot is increased, tonic spiking solutions arise (Fig. 6B). A further increase in c tot accelerates the tonic spiking in v (Fig. 6C), and both ca i and na i oscillate at higher values than in (A) and (B). A graphical summary of the effect of c tot variations on the trajectories is provided in Fig. 7, where the bifurcation structure of the fast-slow subsystem with respect to c tot for l = 0.94 is displayed. We plot the standard L 2 norm of the solution against c tot . The fast-slow system exhibits tristability between bursting (red solid) and two tonic spiking solutions (blue solid) for small c tot values. As c tot is increased, first the bursting branch becomes unstable at a saddle-node bifurcation and then one spiking branch does the same. At the start of an SB cycle, the trajectory is attracted by the stable bursting branch with c tot ≈ 0.6 (red diamond). Solution behavior switches to tonic spiking if c tot is increased. Note that the lower L 2 norm corresponds to larger ca i and na i values; hence, as the trajectory gets attracted by the lower branch of spiking after the saddle-node bifurcation of the upper one, the tonic spiking solution occurs at larger ca i , na i values.
Although Fig. 7 is suggestive, it remains to study more carefully the slow and superslow dynamics in order to understand the mechanisms underlying the switch from bursting to tonic spiking. To do so, we will average over the fast subsystem oscillations. For convenience, we refer to the two regions of (na i , ca i ) space separated by the curve of HC bifurcations that terminates each regular burst as the silent and active regions, respectively (Fig. 5C).
During each interburst interval within the regular bursting epoch, the full model dynamically collapses to a lower-order system governed by the slow variables ca i and na i restricted to S, the manifold of equilibrium points of the fast subsystem. Each interburst interval occurs when the trajectory projected to (na i , ca i ) space lies in the silent region. During the spiking phase of each regular burst, the solution trajectory is still largely determined by the slow variables ca i and na i , but these variables are perturbed by the voltage spike and the Ca 2+ -influx associated with each action potential. This spiking phase corresponds to the active region of (na i , ca i ) space. In this region we employ the method of averaging by numerically averaging the derivatives of the slow variables over one cycle of the action potential, while the superslow variables c tot and l are treated as static parameters. By doing so, we reduce the fast-slow subsystem to two equations for just the slow variables. For g 1 and g 2 defined as the right-hand sides of (5d) and (5f), respectively, the reduced system can be written as We refer to the reduced problem (6a)-(6b) as the averaged slow system. The nullclines of the averaged slow system are curves of (na i , ca i ) values along which there exist periodic solutions (with period T (ca i , na i )) of the fast-slow subsystem that satisfy the additional constraint of either ċ a i = 0 or ṅ a i = 0. In future discussions of the dynamics of the averaged slow system, we will refer to the ca i and na i average nullclines as ca av and na av , respectively. Each intersection of ca av and na av is a fixed point of system (6a)-(6b) representing a tonic spiking solution of the fast-slow subsystem, which we will refer to as FP avi for some index i. Figure 8 illustrates phase planes of the average slow system (6a)-(6b) for l = 0.94 and c tot = 0.6, 0.7, 0.74 as in Fig. 6. In each panel of Fig. 8, the green curve represents the HC bifurcation of the fast subsystem that forms the boundary of the oscillation region. Above HC, where the fast subsystem oscillates (Fig. 5C), the averaged nullclines ca av (blue curve) and na av (green curve) are shown. As noted before, fixed points of (6a)-(6b), FP avi (yellow diamonds), are given by the intersections of these nullclines, and one can usually determine the stability of the fixed points by considering the nullcline configuration.
In Fig. 8A with c tot = 0.6, the two average nullclines intersect at a stable fixed point FP av1 (yellow diamond), which corresponds to the upper spiking branch in Fig. 6D. Despite the existence of this stable fixed point (corresponding to stable tonic spiking), the fast-slow subsystem exhibits bursting since our chosen initial values lie in the basin of attraction of the bursting branch. Correspondingly, in Fig. 8A, the projected trajectory moves clockwise, exhibiting small loops corresponding to spikes within a regular burst, until it crosses HC, at which point the regular burst terminates and the loops are lost while the trajectory transits along a stable branch of the equilibrium curve S (not shown here).
At c tot = 0.7, the stable bursting branch has been lost (Fig. 7) and hence the trajectory is now attracted by the stable fixed point FP av1 (Fig. 8B, yellow diamond). There are also a saddle equilibrium FP av2 , visible in the figure, and a third fixed point of (6a)-(6b) that lies at larger ca i and na i values, not shown here. As a result, the fast-slow subsystem converges to the lower stable fixed point FP av1 and exhibits tonic spiking.
As c tot increases further to 0.74, the lower two fixed points FP av1 and FP av2 collide and annihilate through a saddle-node bifurcation (SN 2 in Fig. 7; note their absence in Fig. 8C) and only the upper stable fixed point FP av3 remains (Fig. 8D, yellow diamond), corresponding to the lower spiking branch in Fig. 7. Therefore, the trajectory  Fig. 8D). Once there, the trajectory approaches the fixed point FP av3 since ṅ a i remains positive. As the trajectory converges toward FP av3 , tonic spiking dynamics at large ca i and na i values results.
Using slow averaged dynamics, we have elucidated how the transition from bursting to tonic spiking occurs and why ca i jumps to larger values as c tot increases for l fixed. Similarly, we summarize the effects of c tot on the average slow system by using bifurcation analysis as shown in Fig. 9. Again, the upper (resp. lower) branch of the bifurcation diagram in L 2 norm corresponds to lower (resp. upper) values of ca i as well as na i and hence solutions on this branch denote FP av1 (resp. FP av3 ). The middle branch represents the unstable saddle FP av2 (see Fig. 8B). Notice that Fig. 9 looks qualitatively the same as the tonic spiking curves shown in blue in Fig. 7, hence either can be used to illustrate the influence of c tot on the oscillatory trajectories for a fixed value of l. As c tot is increased, calcium jumps up at a saddle-node (SN) seen both in the tonic spiking branch in Fig. 7 and in the curve of average system fixed points in  Fig. 9. On the other hand, the onset of spiking happens at a SN of the bursting branch (Fig. 7, red). Next we extend this bifurcation analysis and examine the dependence of the solution patterns of the fast-slow subsystem on both superslow variables c tot and l. To do this, we compute two-parameter bifurcation diagrams in (c tot , l)-space.
In Fig. 10C, the fast-slow subsystem spiking/bursting boundary (solid red, SN 1 ) was calculated by following in the two superslow variables (c tot , l) the SN point where the bursting branch loses stability (Fig. 7, red curve). Also shown is the boundary (solid blue, SN 2 ) demarcating where ca i jumps up, computed by following the upper fold point of the tonic spiking branch in Fig. 7. We can also use direct simulation of the fast-slow subsystem (e.g., fixing c tot , varying l systematically, and then repeating for a different c tot ), to estimate the bursting/spiking boundary curve and the onset of the jump-up in calcium (shown in dashed red and blue). The dashed curves approach the solid ones as we exaggerate the timescale separation between the fast-slow subsystem and the superslow variables (data not shown).
Combining these results with the fast subsystem bifurcation analysis as illustrated in Fig. 5C, and reproduced in Fig. 10B, we are now able to fully understand the SB solution shown in Fig. 10A. Starting from the yellow star at low c tot and l (Fig. 10), a sequence of small bursts is produced as the trajectory in (ca i , na i )-space oscillates between the LF and HC curves multiple times (Fig. 10B). A key step toward termination of this process occurs when the increases of c tot and l push the trajectory in (c tot , l)-space across the SN 1 , such that the fast-slow subsystem enters the tonic spiking regime. In fact, several more small bursts actually occur after the crossing of SN 1 and are followed by the start of tonic spiking at the triangle. In the singular limit, however, this additional bursting will be lost and the tonic spiking occurs when the trajectory reaches SN 1 . Within the same projection, the trajectory evolves rightward from the triangle and eventually passes the SN 2 curve at some point close to the yellow circle, which initiates the jump-up of ca i . The increase in ca i as well as in na i from the yellow circle, which occur on the slow timescale, can be seen in the projection shown in Fig. 10B. Note that this jump corresponds in the projection into (ca i , c tot , l)-space to the convergence of the trajectory to the right branches of the ca i -nullsurface family (Fig. 4B). For ca i large, the trajectory in (ca i , c tot , l)space lies above the l-nullsurface (not shown here); consequently, l next decreases (corresponding also to the decrease in l in Fig. 10C), leading to the reduction of ca i . Therefore, the trajectory in (na i , ca i )-space falls back from the peak in ca i direction  Fig. 7, denoted as SN 1 and SN 2 , respectively. We use SN 1 (solid red) to approximate the bursting/spiking boundary for the fast-slow subsystem (dashed red) and use SN 2 (solid blue) to approximate the onset of jumping up of ca i (dashed blue), respectively and moves towards HC (Fig. 10B). Once it crosses the HC bifurcation curve at the yellow square, the long burst ends and the trajectory enters the silent phase. As the solution finally returns to its starting point (yellow star), one cycle of the SB solution is completed. Remark 1 Some time after ca i jumps up at the yellow circle, the amplitude of the v spikes exhibits a sudden decrease followed by a gradual increase (Fig. 10A). This behavior arises because periodic orbits of the voltage compartment in the Jasinski model initiate in an AH bifurcation with zero amplitude, while orbit amplitudes increase closer to the HC bifurcation (Fig. 4A). Therefore, the sudden decrease of the amplitude of v spikes results from the fact that the jump-up of ca i pushes the trajectory away from the HC curve and closer to the AH. The subsequent decrease in ca i yields a return toward the HC, leading to the gradual increase in the amplitude. This mechanism is the same as observed in MB solutions in [15].

Identifying Timescales
MB solutions, studied previously [15], involve gradual transitions between two different types of bursts, like the SB solutions that we are now considering but with different underlying biological mechanisms. Surprisingly, we obtained the non-intuitive result that the existence of robust MB solutions does not require a third timescale. Thus, a natural question is: how many timescales are fundamentally important for generating SB solutions? To address this question, we adopt the approach used in [23] of transforming our original system into certain two-timescale systems by adjusting system parameters (see Table 5 in Appendix 1). Then we consider whether SB solutions can persist under these adjustments.
The Jasinski model has 7 fast, 2 slow, and 2 superslow (7F, 2S, 2SS) variables. In theory, the timescale separation between some of these groups might not be necessary to generate SB dynamics. Thus, we will consider what happens if we group together the fast and slow variables to form a (9F, 2SS) system and what happens if we group together the slow and superslow variables to form a (7F, 4SS) system. To do so, we first choose A = 1, so that l evolves on a comparable timescale to the other superslow variable, c tot . Since parameters that control the timescale for ca i will also affect the timescale for c tot , we introduce a new parameter, β, with default value 1, as a scaling factor specifically for the right-hand side of the ca i Eq. (1d). To change the timescale for na i , we vary α Na . We will form our (9F, 2SS) system by increasing both β and α Na by a factor of 100, and we will form our (7F, 4SS) system by reducing both β and α Na by a factor of 10.
With its original scaling, system (5a)-(5g) generates a SB solution, as shown in Fig. 10A. Increasing A to 1 does not change this solution qualitatively, except that the number of small bursts decreases. That is, as l becomes faster, the trajectory At least one such solution is present for all c tot ; in fact, the stable branches overlap over a small interval in c tot , yielding bistability projected into (c tot , l)-space will reach the spiking/bursting boundary curve (SN 1 ) earlier and hence small bursts give way to a long burst earlier (Fig. 10C). We next compare this SB solution to solutions from the (9F, 2SS) version of the system (5a)-(5g) described above (compare Fig. 10A with Fig. 11A). The fast/slow decomposition method (between (v, y) and (na i , ca i ) for understanding the mechanisms underlying the regular bursts within the SB solution no longer applies to the (9F, 2SS) case, because ca i and na i now evolve on the fast timescale. In contrast to the original (7F, 2S, 2SS) case where the (7F, 2S)-subsystem generates bursting solutions for relatively small c tot and l values (see Figs. 6 and 7), the new 9-dimensional fast subsystem for the (9F, 2SS) case exhibits tonic spiking for all c tot and l values within a complete bursting cycle. The effect of c tot on the fast subsystem trajectories can be determined from Fig. 11B, where the bifurcation diagram of the fast subsystem (9F) with respect to c tot for l = 0.94 is displayed. Two branches of stable solutions are present, both corresponding to tonic spiking, and one or more stable tonic spiking solutions exist for all c tot values, while the bursting branch that was found in the original system no longer exists in this case. The stable branches persist in two-parameter (c tot , l)space for all relevant l (data not shown); therefore, as c tot , l evolve on the superslow timescale, the trajectory remains on these spiking branches and spiking persists for all time, and the system cannot sustain a SB solution.
Under the alternative rescaling to a (7F, 4SS) system, one cycle of the bursting solution is as shown in Fig. 12A. Only long sigh-like bursts now occur, without regular bursts. Note that the (7F, 2S, 2SS) and (7F, 4SS) systems have the same fast subsystem, with the same LF and HC curves in (na i , ca i )-space (Figs. 10, 12B). Within the projection into the (na i , ca i )-space, a burst of activity begins as the trajectory evolves clockwise in the direction of increasing ca i and na i from the LF curve in the lower left part of the figure. Eventually, c tot , l change enough to cause a rise in the target value of ca i ; the details differ from the (7F, 2S, 2SS) case because the slow averaged dynamics and fast-slow subsystem are no longer relevant, but the outcome is similar. With the (7F, 4SS) rescaling, ca i , na i evolve on the same superslow time scale as c tot , l. Hence, the drift of the trajectory before this transition is too slow for the solution to reach the curve of HC bifurcations and fall silent. As a result, the single burst continues all the way up until the transition; that is, regular bursting never occurs. While the burst continues, all four superslow variables increase until the trajectory projected to (ca i , l, c tot )-space goes above the l-nullsurface (not shown here). After that, l starts decreasing, which eventually leads to the decrease in ca i (Fig. 12B,C). Again similarly as before, this reduction in ca i is fundamental in terminating the burst as it brings the trajectory across the curve of HC bifurcations (Fig. 12B). Afterwards, the solution enters the silent phase and goes back to the starting point, completing one period consisting simply of one single long burst.
In summary, neither of these two-timescale systems, despite the fact that they are the ones with the most similarity to the full three-timescale system, captures the full features of the SB solution shown in Fig. 10A. It is critical that ca i , na i are distinctly slower than the fast voltage and other variables and faster than c tot , l for regular bursts between sighs to occur. We thus conclude that presence of three timescales is necessary for the emergence of the type of the SB solution we have studied, which differs from what was obtained in the pre-BötC MB model [15].

Sigh-Like Spiking in a Single Pre-BötC Neuron
In Sect. 3, we have used the fast-slow decomposition method, bifurcation analysis and the averaging method to understand the mechanisms underlying the SB solution pattern, from which we discover that the short, eupneic bursts are of the well-known square-wave or fold-homoclinic type [16]. A square-wave burster is a kind of relaxation oscillator except that its active state comprises a fast oscillation rather than a quasi-steady plateau. For a square-wave burster, we can convert the upper active or spiking state to a quasi-steady state by removing specifically the fast-spiking components from the model. Following this logic, since the SS solution is composed of a large spike emerging periodically from a pattern of ongoing regular spikes, we may think of the SS solution as a simplification of the SB solution in which fast-spiking components are removed from all active states.
In fact, the discussion in Sect. 2.2 suggests that instead of working with the full SS pattern containing eupneic spikes ( Fig. 2A), we can focus entirely on the activity of the sigh compartment (Fig. 2B). Thus, we simply need to analyze what happens during the prolonged silent phase and how the transition to the large, sigh-like spike occurs, paralleling the analysis in Sect. 3. We will first analyze what timescales are present in the model in the regime that supports the SS solutions of Fig. 2B and how they should be grouped to apply geometric singular perturbation theory. After applying these methods, we will investigate how many timescales are truly required in (3a)-(3e) in order to obtain these SS solutions.
We begin with nondimensionalization, defining new dimensionless variables (v, ca i , c tot , τ ) and voltage, calcium and time scales Q v , Q c and Q t : note that h and l are already dimensionless in (3a)-(3e). Details of the nondimensionalization procedure, including the determination of appropriate values for Q v , Q c and Q t , are given in Appendix 2. From this process, we obtain a dimensionless system of the form wheref 1 ,ĥ 1 ,ĝ 1 ,ĝ 2 ,f 2 ,ĝ 3 are O(1) functions and ε, δ are small parameters. d has size O(1) for c small (during the silent phase) and approximately 1/δ for c large (during the active phase). From this nondimensionalization result, we conclude that v and ca i evolve on a fast timescale, and h and l evolve on a slow timescale. c tot , however, has different evolution rates at different times: it evolves on a slow timescale for c large becauseĝ ≈ δĝ 1 +ĝ 2 is an O(1) function, and it evolves on a superslow timescale for c small sinceĝ ≈ δ(ĝ 1 +ĝ 2 ) is approximately an O(δ) function.
The SS system features bidirectional coupling between v and ca i , which is similar to that in the Jasinski model; therefore, we may also be able to explain the mechanisms underlying SS dynamics from the perspective of how the (v, h) system is driven by the (ca i , c tot , l) system in analogy to our approach for SB solutions in Sect. 3. In the absence of coupling, the (v, h) system has a unique stable equilibrium and no stable oscillatory solution (see Fig. 13A). When the coupling is restored by making g Ca and g CAN nonzero as given in Table 4, the increase of ca i moves the v-nullcline to the upper left (Fig. 13B, red). Meanwhile, the two folds of the nullcline meet and disappear as ca i is made larger than roughly 0.0376. As a result, the unique fixed point of the (v, h) system remains stable for ca i → ∞. That is, there is no value at which ca i can be fixed to yield stable oscillations in (v, h). However, although the fixed point in (v, h) remains stable for all ca i , the SS trajectory projected into (v, h)-space jumps away from the family of fixed points, to larger v, after staying nearby for a finite time. Therefore, we see that the onset of v-spike does not result from the variation of ca i pushing the (v, h) system through any bifurcation at which oscillations are born. To understand how the increase of ca i triggers the fast jump in v, we notice that ca i can be considered to be as fast as v, rather than a slow variable, when it is not near its nullsurface and thus slaved to the slower variables l and c tot . This analysis confirms that since ca i jumps up on a fast timescale, the bifurcation diagram with ca i treated as a static parameter as shown in Fig. 13 no longer plays a role. Hence, a different subsystem classification is needed to explain the SS dynamics.
We use GSPT to define an array of subsystems for system (7a)-(7e), following [23]. Introducing a fast time t f = t s /ε and letting ε → 0, we can derive a 2dimensional fast layer problem that describes the dynamics of the fast variables v and ca i , for fixed values of the other variables: We define the critical manifold M s to be the manifold of equilibrium points of the fast layer problem, i.e., Obtaining slow reduced problems is trickier here since c tot has different scalings at different times. During the silent phase where ca i is relatively small, taking singular limits ε, δ → 0 in (7a)-(7e) yields a system that describes the dynamics of the slow variables h and l for fixed values of c tot , with all variables restricted to the surface of M s , subject to the constraintf 1 =f 2 = 0. Following the terminology used in [23], we call this system the slow reduced layer problem and we define the superslow manifold M ss to be the manifold of its equilibrium points, i.e., To describe the dynamics of c tot restricted to M ss , we define a superslow time t ss = δt s and rewrite (7a)-(7e) as a rescaled system with respect to t ss . Taking the singular limits ε, δ → 0 in the rescaled system yields the superslow reduced problem: On the other hand, during the active phase when ca i is relatively large, we havê g(v, ca i ) = O (1). In other words, c tot evolves on a slow timescale. In this case, taking the limits ε, δ → 0 in (7a)-(7e) gives a system that describes the dynamics of all three slow variables h, c tot and l, subject to the constraintf 1 =f 2 = 0. We call system (11a)-(11c) the slow reduced problem.
Finally, we can also define a fast-slow subsystem of fast and slow variables together, as we did in Sect. 3.1. With these timescale groupings and subsystems defined, we can proceed to analyze the mechanisms underlying the SS solutions.

Analysis of the SS Solutions
The GSPT approach starts with a bifurcation analysis of the layer problem (8a)-(8b), which requires a visualization of the set of equilibria of the layer problem, M s , given byf 1 (v, h, ca i ) =f 2 (v, ca i , c tot , l) = 0. Since the phase space for the full nondimensionalized system (7a)-(7e) is 5-dimensional, M s is a 3-dimensional manifold. Sincef 2 does not depend on h, the projection of M s onto (v, h, ca i )space is simply given byf 1 (v, h, ca i ) = 0; that is, for all relevant (v, h, ca i ), we can solvef 2 = 0 for an appropriate choice of (c tot , l). We can solvef 1 = 0 for h as a function of v and ca i and can therefore represent M s projected onto (v, h, ca i )-space as h = F 1 (v, ca i ) for a function F 1 (Fig. 14A). For each fixed ca i value,f 1 = 0 is represented by a single curve, as shown in Fig. 13B.

Remark 2
For the purpose of nondimensionalization in Appendix 2, we require that c tot ∈ [0, 1]. However, technically, the value of c tot needed to makef 2 (v, ca i , c tot , l) = 0 for some relevant (v, h, ca i ) can be greater than 1, but only as large as 1.5. Since this number is just marginally larger than 1, adjusting our nondimensionalization in Appendix 2 accordingly would not change the scales of any variables, and hence this is just a technicality that can be ignored for practical purposes.
Further insight into the dynamics of calcium comes from viewing the trajectory and M s in (ca i , c tot , l)-space. To visualize this projection, we notice that during the long silent phase of the SS solution, v restricted to the surfacef 1 = 0 is approximately a constant denoted byv ≈ −0.5 (see Figs. 2C, 13B and 14A).f 2 does not depend on h, so we can substitutev ≈ −0.5 and solvef 2 (v, c, c tot , l) = 0 for c tot as a function of c and l and can therefore readily visualize the projection of M s onto (c, c tot , l)space (Fig. 14B). Also shown in Fig. 14A and B are the projections of the solution trajectory (black) and M ss (red). As v jumps up to a larger range of values during the v-spike, the critical manifold M s should no longer be represented in (ca i , c tot , l)space by the ca i -nullsurface fixed atv. However, reasons similar to those discussed in Sect. 3, the calcium equation only depends weakly on v during the active phase; that is, as the trajectory enters the spiking phase, the ca i -nullsurfaces for v =v and for v = v max lie extremely close to each other. Therefore, it suffices to consider the single ca i -nullsurface for v = −0.5 to qualitatively understand the projection of the full SS dynamics shown in Fig. 14B.
By having visualizations of M s and M ss projected onto both (v, h, ca i )-space and (ca i , c tot , l)-space, we can now understand the evolution of the SS solution in terms of the shapes and relative positions of M s and M ss . Starting from the yellow star in Fig. 14, the trajectory is in the silent phase and evolves on the slow timescale under the slow reduced layer problem (9a)-(9b), until it approaches sufficiently close to the superslow manifold M ss (near the yellow circle). From there, the trajectory evolves on the superslow timescale under (10). It follows M ss to a subcritical AH bifurcation (yellow diamond) of the fast-slow subsystem with respect to bifurcation parameter c tot , where a branch of unstable small amplitude periodic orbits is born (not shown here) and M ss becomes unstable (yellow diamond). From there, the trajectory makes a fast jump to large v and ca i , governed by (8a)-(8b); this jump corresponds to the onset of a spike in v in the SS solution. After the trajectory reaches M s after the fast Fig. 15 Simulation of the sigh-like spiking solution. Sigh-like spiking solution of (7a)-(7e). The yellow symbol indicates the transition point between slow and superslow flow jump, the active phase begins. As we can see from Fig. 14B, ca i is relatively large near the right branch of the ca i -nullsurface and therefore, as discussed previously, c tot is no longer a superslow variable. Thus, during the spiking phase, the flow is governed by the slow reduced problem (11a)-(11c). Since there are no branches of M ss on M s for ca i large, the trajectory moves on the slow timescale along M s under (11a)-(11c) until it meets the lower fold of the ca i -nullsurface (Fig. 14B), from which it jumps back to its starting point (yellow star) on the fast timescale under (8a)-(8b), completing a full cycle.
Based on the above analysis, we summarize the different timescales on which v evolves in Fig. 15, where the yellow circle indicates the transition point between the slow and superslow timescales that occurs as the trajectory reaches a small neighborhood of M ss as it evolves along the lower-v surface of M s .

Identifying Timescales
Next, we seek to identify whether the SS solution discussed above is truly a threetimescale phenomenon. In our original scaling, the Toporikova model has 2 fast, 2 slow and 1 superslow (2F, 2S, 1SS) variables. Similarly to Sect. 3.2, we assess the importance of having three timescales in two natural ways, by adjusting the two slow variables to be either fast or superslow. That is, we first speed up h and l by decreasinḡ τ h and increasing A by a factor of 100, respectively, so that they evolve on the same timescale as v and ca i , resulting in a 4 fast, 1 superslow (4F, 1SS) system. Second, we slow down h and l by adjustingτ h and A in the opposite way to produce a system with 2 fast and 3 superslow (2F, 3SS) variables. Then we consider whether or not these two-timescale systems can generate solutions that are similar to the SS solution.
In the (4F, 1SS) rescaling, a different type of trajectory lacking large v and ca i spikes is observed (see Fig. 16A and B for different projections of the solution). For the fast layer dynamics of the (4F, 1SS) system, the critical manifold is our former superslow manifold, M ss , which lies within M s . The outer branches of M ss are stable while the middle branch is unstable with respect to the layer system. The stable trajectory of the (4F, 1SS) system is attracted to a stable branch of M ss . Eventually, the trajectory passes the fold of M ss (blue triangle) where M ss destabilizes, a transition to the fast layer problem occurs, and hence the solution of the (4F, 1SS) system follows the fast layer flow to another stable branch of M ss . The subsequent motion is governed by the superslow flow until v and ca i jump down from another fold of Under the alternative rescaling to a (2F, 3SS) system, the projections of the solution trajectory are as shown in Fig. 16C and D. This system has the same critical manifold M s as the original system. In contrast to the (4F, 1SS) system, M ss , although shown for comparison, is no longer meaningful in the (2F, 3SS) system. Hence, the trajectory simply follows M s , and a standard two-timescale relaxation oscillation results.
In summary, neither of these two-timescale systems captures the full features of the oscillating solution shown in Fig. 2C. Since a two-timescale system will lose either the large spikes in the v and ca i time series or the two-scale aspect of the recovery of v between spikes, we classify the SS solution as an intrinsically three-timescale phenomenon. We note, however, that the qualitative difference between the SS solution for the (2F, 2S, 1SS) and (2F, 3SS) scalings are not nearly as significant as the differences arising in the three-and two-timescale scalings for the SB model (1a)-(1g).

Comparison of Results for SB and SS Solutions
We are now in a position to compare mechanisms underlying SB and SS solutions and how various timescales are manifested in the two solutions. At a general initial condition within the equilibrium set of the relevant fast subsystem (denoted by S in the Jasinski model and M s in the Toporikova model, respectively), both SB and SS solutions evolve on a slow timescale. From such an initial condition, the solutions will converge to different dynamic regimes. Specifically, the SB solution is attracted by a stable bursting state for the fast-slow subsystem (as shown in Fig. 7), yielding multiple regular bursting cycles. In contrast, the SS solution converges to M ss and hence enters the quiescent state, where a plateau of v results. Although constrained to different fast-slow subsystem attractors, both solution trajectories are governed by the dynamics of the superslow variables, variations of which will then push the fast-slow subsystem, which consists of all the fast and slow variables with superslow variables as static parameters, through some bifurcations at which either a sigh-like burst or a sigh-like spike can occur.
Despite this commonality, the mechanisms underlying the transitions from regular bursts to the sigh-like burst and from the quiescent plateau to the sigh-like spike differ in the two models. Specifically, for the SB solution, the sigh-like burst arises as the increase of superslow variables, c tot and l, switches the behavior of the fastslow subsystem from bursting to tonic spiking as the bursting branch loses stability at a SN bifurcation of periodic orbits (SN 1 in Figs. 7 and 10C). In contrast, the sighlike spike in the SS solution occurs as the trajectory crosses an AH bifurcation of the fast-slow subsystem of (7a)-(7e), where M ss becomes unstable. Different bifurcations associated with the generation of sighs in the two models induce different requirements on the timescale of calcium to produce sighing activity. As discussed in Sect. 3.2, ca i cannot be grouped together with v as a fast variable to generate SB solutions. There is no such constraint on ca i in producing SS solutions, however. If an experimental manipulation could accelerate calcium dynamics, then the loss of sighing dynamics would support the SB model, which involves a unified mechanism for eupneic and sigh-like burst generation, whereas little change in sighing dynamics would support the SS model, which is based on separate eupneic and sigh-like burst generation mechanisms.
Another difference between SB and SS solutions lies in the transition of ca i to large values. In the SS solution, ca i jumps up after the AH bifurcation on M ss , which is the same bifurcation inducing the sigh-like spike. The jump-up of ca i in the SB solution, however, happens at a different bifurcation than that for the onset of a sighlike burst. After an additional drift on the superslow timescale, the transition occurs at a SN bifurcation of a spiking branch of the fast-slow subsystem of the Jasinski model (SN 2 in Figs. 7 and 10). In other words, the observation that a surge in calcium coincides with the onset of sighing would support the SS model, whereas a finding that the surge in calcium occurs after sighing is already under way would support the SB model.
After the rapid increase of ca i , the trajectory enters the active phase, during which the SB model generates continuous spiking constituting the sigh-like burst, which is still governed by the dynamics of superslow variables (i.e., c tot and l). In the SS model, during the active phase when the big spike occurs, the evolution of the solution along the right branch of the calcium nullsurface is governed by the slow reduced problem and hence occurs on a slow timescale, such that the SS sigh event has a shorter duration than the sigh event in the SB solution.
Finally, termination mechanisms for both SB and SS solutions involve ca idependent inactivation of [IP 3 ] (regulated by l), resulting in the reduction of ca i and deactivation of I CAN [13]. In both SB and SS solutions, regular oscillations (bursting/spiking) persist when I CAN is removed, whereas the sigh-like oscillations are eliminated under the supression of I CAN or I Ca , confirming that the sigh-like activities in both models critically depend on both I CAN and I Ca .

Discussion
To understand the mechanisms underlying the generation of sighs, we considered two distinct single-compartment models for respiratory neurons in the pre-BötC that have the ability to generate sigh-like activity. For both models, we performed nondimensionalization to identify possible groupings of variables into classes corresponding to they timescales on which they evolve. After establishing such groupings, in some cases with class membership changing depending on the location of the trajectory in phase space, we used a non-rigorous GSPT approach to elucidate the roles of various timescale-based subsystems and their bifurcation structures in producing sigh-like dynamics. We identified commonalities and differences between the mechanisms involved for the two models and argued that for both models, the sigh-like dynamics involves three timescales in an essential way. This work adds to the growing literature of dynamical systems analyses of three-timescale systems [15,23,[25][26][27][28][29][30] and in particular to recent efforts to identify how many timescales are really needed to produce particular dynamic patterns [15,23]. It also provides information about model parameter relationships needed to support sigh-like activity, which may be useful for future efforts to model the repertoire of pre-BötC dynamics and their variations under normal conditions, environmental and metabolic challenges, and pathologies [2].
The first model that we considered is a self-coupled neuron model featuring I NaP , I CAN and the Na + /K + -pump current. An aspect of this model that is more biologically realistic than previous models studied in [15] and [23] is the inclusion of bidirectional coupling between voltage and calcium dynamics. In Sect. 3, we extended and applied analysis methods from [15] to the Jasinski model (1a)-(1g) and explained the mechanisms underlying its SB solutions. While the bidirectional coupling between V and Ca i as well as more detailed Ca 2+ dynamics make the implementation of the decomposition method more difficult than in past work on similar models, fast-slow averaging allowed us to complete the analysis. Besides describing specific details of the SB solution features, we have also investigated whether this solution fundamentally involves three timescales. Unlike the MB solution in [13][14][15], our analysis shows that SB solution features are lost under the natural groupings to two timescales, supporting the preliminary conclusion that SB dynamics in system (1a)-(1g) requires at least three timescales. A more rigorous demonstration of this requirement is still an open matter, and indeed rigorous proofs that particular solution types can only occur when three (or more) timescales are present have not, to our knowledge, been provided in the literature to date.
Several conditions that support the existence of the SB solutions can also be deduced from our analysis. To obtain SB solutions, we require relatively small g Ca , whereas large g Ca will eliminate SB patterns in the Jasinski model. That is, the increase of g Ca will speed up Ca i and Ca tot (see Table 5) such that the time available for the regular bursting phase becomes shorter and as a result, the number of small bursts decreases. With further increases in g Ca , the regular bursts will completely disappear and the SB solution will be lost. Second, the regular bursting phase also necessitates a long enough time before Ca i jumps up to allow the solution trajectory to undergo multiple crossings between the LF and HC curves in the slow (Na i , Ca i )-space. As suggested by the (9F, 2SS) case shown in Fig. 11, each crossing between the LF and HC curves should also be slow enough for the full system to generate a burst, rather than simply a spike. On the other hand, if the evolution of Na i and Ca i are too slow for the solution to complete a single square-wave burst before the fast subsystem transitions to tonic spiking due to the evolution of Ca tot and l, as shown in Fig. 12, the SB solution no longer exists. This suggests that the existence of the SB solution requires the timescale for Ca i , Na i to be faster than Ca tot , l. These types of arguments could be useful for deriving a minimal biological model for SB dynamics, although we do not complete this step here since we do not have a particular application in mind for such a model. Similar analysis can also been used to adjust parameters in order to enhance the robustness of SB solutions with respect to other parameter variations, as has been done for another multiple-timescale respiratory neuron model [15].
Another possible direction for future work is to explain why the self-excitation in system (1a)-(1g) appears to be required for the SB solution to exist. In our setting, we have chosen g SynE = 20 to obtain the SB solutions and a sufficient decrease of g SynE will eliminate the SB dynamics, which agrees with the findings from [8] that SB behavior requires excitatory synaptic inputs. A preliminary numerical simulation shows that a decrease of g SynE can result in a transition from bursting to spiking in the voltage compartment, preventing the regular square-wave bursting phase from occurring at low Ca i as needed for SB dynamics.
In addition to the Jasinski model, we have also considered a second model for inspiratory pre-BötC neurons (3a)-(3e), which can yield SS solutions [9]. Similarly as for the Jasinski model, we applied GPST to understand the dynamics of the SS solution. By doing so, we discovered that the SS solution is not just a simple reduction of the SB solution in the same way that a relaxation oscillator can be viewed as a simplification of a square-wave burster. Instead, the models differ in terms of the mechanisms underlying their high-amplitude sigh-like activities as well as in the details of which timescales control particular solution features, as discussed in Sect. 5. Nonetheless, sigh-like activities in both models depend critically on calcium oscillations, consistent with the previous experimental data showing that I Ca blocker and I CAN blocker [6,7,24] could terminate sighs generated in medullary slices containing the pre-BötC without suppressing regular bursting activity.
As with the SB solution, we also demonstrate that SS dynamics appears to be a three-timescale behavior. In addition, an interesting mathematical problem arising from the analysis of the SS solution is the development of systematic methods to treat equations that evolve on different timescales in different parts of phase space (cf. [31]). In Sect. 4, we have illustrated in a non-rigorous way how the GPST approach can be extended to a multiple-timescale model where the timescale of one variable, Ca tot , changes with respect to the location of a trajectory in phase space.
Moreover, according to the nondimensionalization results in Appendix 2, we find that the calcium derivative, denoted by f 2 , can also be represented as a combination of functions with different relative sizes. Therefore, in different parts of phase space, Ca i can evolve on the fast, slow or superslow timescale. We repeated our GSPT analysis in a way that takes this observation into consideration, but no qualitatively additional information was gained from doing so. Hence, Ca i is simply treated as a fast variable in our analysis of SS dynamics in this paper. Development of a more systematic and rigorous approach for assessing when phase-dependent scalings are important would be a helpful step for future analyses.
The dimensionless currentsĪ x appearing in (12a) and (12f) are given by I x g max Q v . Dimensionless calcium fluxes are given byJ ER IN = J ER IN σ/P max andJ ER OUT = J ER OUT σ/P max , where P max takes the value 1,148 according to numerical results. In Eqs. (12c) and (12d), K c tot = τ Ca · α Ca · g Ca · Q v /Q c and K c = K Ca ·P max α Ca ·σ ·g Ca ·Q v . Since we expect V ∈ [−60, 20], a suitable choice for the voltage scale is Q v = 100 mV. Similarly, we choose Q na = 30 mM since Na i ∈ [16,21]. Q c , the scale for Ca i and Ca tot , should be determined by magnitudes of both variables. Numerical simulations show that Ca i varies from O(10 −5 ) to O(10 −3 ), while Ca tot is roughly O(10 −3 ) ( Fig. 4B and C). Therefore, we choose Q c = 10 −3 mM. Using these scales, we see that all terms in the right-hand sides of Eqs. (12a)-(12g) are bounded (in absolute value) by one except that the last term in (12d) is roughly O (10). To resolve this issue, we divide both sides of (12d) by 10 and obtain the following: the right-hand side of which is now O(1), as we desire for our nondimensionalization. The coefficients of the derivatives in the left-hand sides of Eqs. (12a)-(12g) with (12d) replaced by (13) now reveal the relative speeds of evolution of the variables. We summarize the timescales for variables of the Jasinski model in Table 5.
Comparing the values of these coefficients indicates how fast each corresponding variable is; the larger the value, the slower the corresponding variable. According to our nondimensionalization results summarized in Table 5, we choose to group together V , the gating variables m Na , h Na , m Ca , h Ca , m K , and s as the variables evolving on a fast timescale. We group (Ca i , Na i ) as evolving on a slow timescale and (Ca tot , l) as evolving on a superslow timescale. We choose the slow timescale as our reference time, i.e., pick Q t = 100 ms and let R x denote coefficients of dx dt in equations (12a)-(12b), (13) and (12e)-(12g), the dimensionless system then becomes the system (5a)-(5g) given in Sect. 3.1, namely, with dimensionless currentsĪ x = I x g max Q v . We have now nondimensionalized the equations for V and h, and next we deal with equations for the calcium compartment (15c)-(15e). From numerical simulations, Ca i and Ca tot typically lie between 0 μM and 1 μM. Therefore, we define G c = max(G 3 (Ca i )), G S = max(g SERCA (Ca i )) and Φ = max(φ(Ca i )) over the range Ca i ∈ [0, 1] and define P max to be the maximum of {L IP 3 , P IP 3 G c , G S }. From system (15a)-(15e), we obtain the following dimensionless system: Since we expect V ∈ [−60, −20] and Ca tot , Ca i ∈ [0, 1], suitable choices for the voltage and calcium scales are Q v = 100 mV and Q c = 1 μM, respectively. We also see that values of h ∞ , h and l all lie in the range [0, 1]. For the choice of parameters specified in Table 4, the maximum conductance is g K , so we have g max = g K . Numerical evaluation of 1/τ h (V ) for V ∈ [−60, −20] shows that T h ≈ 8 × 10 −4 ms −1 = O(10 −3 ) ms −1 . We also obtain G c ≈ 0.0456 and G S ≈ 1,000 pL · ms −1 , so we have P max ≈ 1412 pL · ms −1 and hence K c ≈ 7 × 10 3 = O(10 4 ). Similarly, we have K c tot ≈ 30. Using these values, we see that all terms on the right-hand sides of (16a)-(16e) are bounded (in absolute value) by one except K c tot and K c . To resolve these issues, we divide both sides of (16c) and (16d) by K c tot and K c , respectively, and keep other equations unchanged from (16a)-(16e):