Krasimira Tsaneva-Atanasova
email: k.tsaneva-atanasova@bristol.ac.uk

Hinke M Osinga

Thorsten Rieß

Arthur Sherman

Full system bifurcation analysis of endocrine bursting models

Keywords: Excitable systems, Bifurcation theory, Bursting oscillations, Spike Adding, Endocrine Cells

Plateau bursting is typical of many electrically excitable cells, such as endocrine cells that secrete hormones and some types of neurons that secrete neurotransmitters. Although in many of these cell types the bursting patterns are regulated by the interplay between voltage-gated calcium channels and calcium-sensitive potassium channels, they can be very different. We investigate so-called square-wave and pseudo-plateau bursting patterns found in endocrine cell models that are characterized by a super-or subcritical Hopf bifurcation in the fast subsystem, respectively. By using the polynomial model of Hindmarsh and Rose (Proceedings of the Royal Society of London B 221(1222), 87-102), which preserves the main properties of the biophysical class of models that we consider, we perform a detailed bifurcation analysis of the full fast-slow system for both bursting patterns. We find that both cases lead to the same possibility of two routes to bursting, that is, the criticality of the Hopf bifurcation is not relevant for characterizing the route to bursting. The actual route depends on the relative location of the full-system's fixed point with respect to a homoclinic bifurcation of the fast subsystem. Our full-system bifurcation analysis reveals properties of endocrine bursting that are not captured by the standard fast-slow analysis.

Introduction

Plateau bursting is an intrinsic property of endocrine cells by means of which they achieve the rise in intracellular calcium concentration necessary for secretion [START_REF] Barg | Insulin secretion: a high-affinity Ca 2+ sensor after all?[END_REF][START_REF] Beauvois | Glucose-induced mixed [Ca 2+ ] c oscillations in mouse betacells are controlled by the membrane potential and the SERCA3 Ca 2+ -ATPase of the endoplasmic reticulum[END_REF][START_REF] Jonas | Temporal and quantitative correlations between insulin secretion and stably elevated or oscillatory cytoplasmic Ca 2+ in mouse pancreatic 𝛽-cells[END_REF][START_REF] Van Goor | Dependence of pituitary hormone secretion on the pattern of spontaneous voltage-gated calcium influx. cell type-specific action potential secretion coupling[END_REF]. It is mediated via the interaction of voltage-gated calcium (Ca 2+ ) channels and various potassium (K + ) channels in the cell membrane. Upon stimulation that generally leads to depolarization, the membrane potential becomes more positive, opening the voltage-gated Ca 2+ channels. The resulting Ca 2+ influx into the cytosol triggers activation of calcium-sensitive potassium (K Ca ) channels, generating the outflow of K + ions that repolarizes the membrane potential. This repolarization leads to closure of voltage-gated Ca 2+ channels and subsequent decrease in the cytosolic calcium levels ([Ca 2+ ] i ). This sequence of events leads to oscillations in the intracellular calcium concentration that are accompanied by plateau-bursting electrical activity. Such prolonged electrical activity is an efficient way to increase intracellular Ca 2+ , in contrast to brief neuronlike single spikes. The increase of [Ca 2+ ] i stimulates the release of hormones from secretory vesicles [START_REF] Barg | Insulin secretion: a high-affinity Ca 2+ sensor after all?[END_REF][START_REF] Beauvois | Glucose-induced mixed [Ca 2+ ] c oscillations in mouse betacells are controlled by the membrane potential and the SERCA3 Ca 2+ -ATPase of the endoplasmic reticulum[END_REF][START_REF] Jonas | Temporal and quantitative correlations between insulin secretion and stably elevated or oscillatory cytoplasmic Ca 2+ in mouse pancreatic 𝛽-cells[END_REF][START_REF] Van Goor | Dependence of pituitary hormone secretion on the pattern of spontaneous voltage-gated calcium influx. cell type-specific action potential secretion coupling[END_REF]. The plateau-bursting electrical activity is characterized by periodic switches between an active (depolarized) phase accompanied by increase in [Ca 2+ ] i and a silent (repolarized) phase during which [Ca 2+ ] i decreases due to Ca 2+ extrusion. Owing to the importance of this activity, numerous modeling studies have been carried out of plateau-bursting in a variety of cell types, including pancreatic 𝛽-cells [START_REF] Bertram | The phantom burster model for pancreatic beta-cells[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Chay | Minimal model for membrane oscillations in the pancreatic beta-cell[END_REF][START_REF] Goforth | Calcium-activated K + channels of mouse beta-cells are controlled by both store and cytoplasmic Ca 2+ : experimental and theoretical studies[END_REF][START_REF] Tsaneva-Atanasova | Diffusion of calcium and metabolites in pancreatic islets: killing oscillations with a pitchfork[END_REF][START_REF] Zhang | The Ca 2+ dynamics of isolated mouse beta-cells and islets: implications for mathematical models[END_REF] and pituitary cells [START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF][START_REF] Van Goor | Paradoxical role of large-conductance calcium-activated K + (BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells[END_REF]. The models in these studies are generally derived using the Hodgkin-Huxley formalism [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] and generate bursting behavior by taking into account the crosstalk between voltage-dependent Ca 2+ channels and K Ca channels in combination with the slow dynamics of intracellular calcium concentration.

An important feature of plateau bursters is that the fast spikes during the active phase do not cross the baseline membrane potential in the silent phase. There are two types of plateau-bursting patterns that have been found in models. The classical square-wave [START_REF] Bertram | Topological and phenomenological classification of bursting oscillations[END_REF][START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF] (or fold-homoclinic (28)) bursting is observed in models for pancreatic 𝛽-cells and is characterized by well-defined spikes in the active phase that correspond to stable limit cycle solutions in the fast subsystem of such models. (Note that 𝛽-cells also exhibit other behaviors that combine with square-wave bursting to produce a variety of oscillations [START_REF] Goel | The geometry of bursting in the dual oscillator model of pancreatic 𝛽-cells[END_REF][START_REF] Pedersen | Phantom bursting is highly sensitive to noise and unlikely to account for slow bursting in 𝛽-cells: Considerations in favor of metabolically driven oscillations[END_REF].) The other type of plateau bursting is typical of pituitary cells [START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF][START_REF] Van Goor | Paradoxical role of large-conductance calcium-activated K + (BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells[END_REF] which exhibit small irregular spikes in the active phase. This pattern has been called pseudoplateau bursting [START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF] because it is produced by transient oscillations rather than stable limit cycles, and has been classified mathematically as fold-subHopf [START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF]. Despite these differences, the bursting patterns observed in endocrine cell models are both governed by the interplay between voltage-dependent Ca 2+ channels and

A c c e p t e d m a n u s c r i p t

K Ca channels. The focus of this paper is investigation of these two plateau-bursting patterns. We are interested in identifying fundamental properties of this class of models in terms of parameter(s) that control the behavior of the system. We consider a generic simplified biophysical model system (see Appendix A.1 for details) based on elements drawn from several published models [START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF]. We show that a model system of the above class can generate both types of plateau bursting, depending on the balance between inward and outward currents. We begin by performing a classical fast-slow bifurcation analysis [START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF] of the generic model and demonstrate how the transition between square-wave and pseudo-plateau bursting takes place in the 'frozen' fast subsystem. We then go beyond the fast-slow analysis to examine bifurcations of the full bursting system with respect to the speed of the slow variable. This is of particular interest for fold-subHopf bursting, which produces appropriate spiking patterns when the separation of time-scales is not extreme.

In order to facilitate the computations and demonstrate the generality of the results, we employ a polynomial reduction of Hindmarsh-Rose type [START_REF] Hindmarsh | A model of neuronal bursting using 3 coupled 1st order differential-equations[END_REF]. We demonstrate that this model, although phenomenological, duplicates the qualitative behavior of the biophysical system in the frozen case and then use it to study systematically the full-system bifurcations that lead to bursting and that control the number of spikes per burst. We follow representative periodic orbits of the full ODE system with different numbers of spikes and find that the same patterns of bifurcations govern both square-wave and pseudo-plateau periodic solutions. Most interestingly, we find for both classes of bursters that the patterns differ depending on the location of the full-system steady state, which lies on a branch of saddle equilibria in the fast subsystem, relative to the homoclinic orbit of the fast subsystem.

Fold-homoclinc bursting has been previously studied and others have similarly pointed out the central role of bifurcations of the full system [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF][START_REF] Wang | Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle[END_REF]. One study [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF] has demonstrated locally the emergence of such bursting oscillations from homoclinic connections in the limit 𝜖 → 0. None of these studies, however, has treated the case of fold-subHopf bursting nor systematically compared it to foldhomoclinic bursting.

Results

Generic Endocrine System

Since the pioneering work of Rinzel [START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF], it has become a standard approach to study bursting oscillations using fast-slow analysis, i.e. by decomposing the

A c c e p t e d m a n u s c r i p t

model into fast and (one-dimensional) slow subsystems and analyzing the dynamics of the full system in the limit of the slow variable treated as a bifurcation parameter. In the class of models (Appendix A.1) that we consider, the only slow variable is the cytosolic calcium concentration (𝑐), which indeed changes much more slowly than the membrane potential (𝑉 m ) and the channel gating variables. The separation of time scales in these models is controlled by the parameter 𝑓 c , which represents the fraction of free [Ca 2+ ] i . Fast-slow analysis has been widely used in theoretical studies of both square-wave and pseudo-plateau bursters [START_REF] Bertram | Topological and phenomenological classification of bursting oscillations[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF], assuming that [Ca 2+ ] i is a slow variable. These studies have shown that square-wave bursting models are characterized by a supercritical Hopf bifurcation in the fast subsystem, whereas in pseudo-plateau bursters this bifurcation is subcritical. We illustrate this in Fig. 1 using the generic endocrine model (Appendix A.1), where we vary the half-maximum activation 𝑉 m L of the voltage-gated Ca 2+ channels in order to change the type of the Hopf bifurcation of the fast subsystem (Eqns. ( 4)-( 5)). We also superimpose the bursting trajectories of the full system for each value of 𝑉 m L on the corresponding bifurcation diagrams.

For 𝑉 m L = -22.5 mV, illustrated in Fig. 1(a), the generic model behaves like a typical model for pancreatic 𝛽-cell square-wave bursting [START_REF] Bertram | Topological and phenomenological classification of bursting oscillations[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF]. Note that the simplified generic model corresponds to classical bursting with cytosolic Ca 2+ as the single slow variable, which is exemplified by the original Chay-Keizer model [START_REF] Chay | Minimal model for membrane oscillations in the pancreatic beta-cell[END_REF]. More recent 𝛽-cell models incorporate effects of Ca 2+ in the endoplasmic reticulum and ATP [START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF], but we neglect these as they have no analog in models of pituitary bursting. As shown in Fig. 1(a), the fast subsystem is bi-stable for a range of values of the control parameter 𝑐, and is characterized by a Z-shaped equilibrium curve that folds at saddle-node bifurcation points labeled SN 1 and SN 2 . The upper branch of this curve consists of stable foci that lose stability as calcium increases via a supercritical Hopf bifurcation (HB). Between HB and SN 1 the upper branch points are unstable foci/nodes that become saddles at SN 1 . These equilibria gain stability again at SN 2 and beyond this point are stable nodes. The Hopf bifurcation gives rise to a branch of stable periodic orbits that terminates in a homoclinic bifurcation (HC) with one of the saddles on the equilibrium branch delimited by the two saddle-node bifurcations SN 1 and SN 2 . Bursting also relies on the fact that the 𝑐-nullcline of the full system Eqns. ( 4)-( 6) (not shown in Fig. 1) intersects the Z-shaped equilibrium curve of the fast subsystem Eqns. ( 4)-( 5) somewhere in the middle branch. Thus, the membrane potential (𝑉 m ) in the full system periodically switches between silent and active phases due to repeated intersections of the bursting trajectory and the 𝑐-nullcline that force it to change direction in the phase space.

For the case 𝑉 m L = -27.5 mV shown in Fig. 1(b), the bifurcation curves

A c c e p t e d m a n u s c r i p t

closely resemble those from several pituitary bursting models [START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF]. Compared with the case for 𝑉 m L = -22.5 mV (Fig. 1(a)), we observe that the 5 mV left shift in the Ca 2+ current activation curve preserves the Z-shaped equilibrium curve of the fast subsystem Eqns. ( 4)-( 5) along with its stability properties, but it has shifted to the right. Furthermore, the Hopf bifurcation of the fast subsystem is now subcritical, which results in a branch of unstable periodic orbits that terminates at HC. As can be seen from the trajectory of the full system Eqns. ( 4)- [START_REF] Belan | Calcium clamp in isolated neurones of the snail Helix pomatia[END_REF] shown in Fig. 1(b), the transition in the type of the Hopf bifurcation results in pseudo-plateau rather than square-wave bursting, as the spikes are due to a slow oscillatory approach to the upper equilibrium branch, not stable limit cycles. Such spikes can only occur if [Ca 2+ ] i is not too slow. Specifically, the rate of increase of [Ca 2+ ] i cannot be much slower than the rate of approach of the solution to the upper equilibrium. Indeed, pseudo-plateau bursters can lose bursting oscillations when the slow variable is made too slow, if the trajectory is absorbed in a stable state on the upper branch of the Z-curve. When the slow variable is faster, however, bursting is possible because the trajectory exits the active phase before reaching the stable steady-state. In the class of plateau-bursting models that we focus on, the parameter 𝑓 c that controls the separation of time scales typically ranges from 10 -3 to 10 -1 [START_REF] Bertram | The phantom burster model for pancreatic beta-cells[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Chay | Minimal model for membrane oscillations in the pancreatic beta-cell[END_REF][START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF][START_REF] Tsaneva-Atanasova | Diffusion of calcium and metabolites in pancreatic islets: killing oscillations with a pitchfork[END_REF][START_REF] Van Goor | Paradoxical role of large-conductance calcium-activated K + (BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells[END_REF], which is only moderately small. For intermediate values of 𝑉 m L (not shown), the bifurcation diagram of the fast subsystem Eqns. ( 4)-( 5) deforms continuously via a (codimension-two) degenerate Hopf bifurcation point, where HB changes from supercritical to subcritical. At first, the branch of unstable periodic orbits turns around at a saddle-node of periodics (SNP) bifurcation, which leads to a branch of stable limit cycles that terminates at HC. As 𝑉 m L decreases and the Z-curve shifts to the right, the point HB also shifts to the right, but to a greater extent and, thus, moves closer to SN 1 and the middle branch of the Z-curve. Hence, eventually the SNP and HC occur simultaneously, after which the periodics never become stable.

The changes in the bifurcation diagrams in Fig. 1 reflect several biophysical effects of the left shift of the Ca 2+ current activation curve, which alters the balance between Ca 2+ and K + currents in the inward direction. This means that more K Ca current is needed to repolarize the bursts, so that the Z-curve shifts to the right. The shift of the HB reflects an enlarged region of conduction block, in which excessive inward current prevents spiking and results in a depolarized plateau. Note that the loss of true spiking increases the Ca 2+ concentration because the mean membrane potential is higher without the hyperpolarized inter-spike interval. The model thus suggests that the levels of intracellular calcium concentration ([Ca 2+ ] i ) during pseudo-plateau bursting could be significantly greater, other things being equal, than during square-wave bursting (Fig. 1

The Polynomial Model

We complement the classical fast-slow analysis with a bifurcation analysis of the full system. Such a full-system analysis provides a different view of the bursting solution as a periodic orbit with a complicated internal structure. This approach is necessary to detect chaos, which, as shown by Terman [START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF], is more robust when the slow variable is not very slow and thus is more likely to be observed in experiments. In order to investigate systematically the full system bifurcation structure of endocrine models we construct a polynomial plateau-bursting model by building into it the common dynamical features found in a number of biophysical modeling studies [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Bertram | The phantom burster model for pancreatic beta-cells[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF][START_REF] Tsaneva-Atanasova | Diffusion of calcium and metabolites in pancreatic islets: killing oscillations with a pitchfork[END_REF].

Equations and Assumptions

The polynomial model is a modified Hindmarsh-Rose type model [START_REF] Hindmarsh | A model of neuronal bursting using 3 coupled 1st order differential-equations[END_REF] with parameters chosen such that the bifurcation diagram of the fast subsystem is similar to that of Eqns. ( 4)- [START_REF] Beauvois | Glucose-induced mixed [Ca 2+ ] c oscillations in mouse betacells are controlled by the membrane potential and the SERCA3 Ca 2+ -ATPase of the endoplasmic reticulum[END_REF], that is, the upper equilibrium branch exhibits a single Hopf bifurcation; compare also [START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF]. The equations have the general form

𝑑𝑥 𝑑𝑡 = 𝑓 (𝑥, 𝑦, 𝑧), (1)
𝑑𝑦 𝑑𝑡 = 𝜙𝑔(𝑥, 𝑦), (2)
𝑑𝑧 𝑑𝑡 = 𝜖ℎ(𝑥, 𝑧), (3)
where 𝑓 (𝑥, 𝑦, 𝑧), 𝑔(𝑥, 𝑦) and ℎ(𝑥, 𝑧) for (𝑥, 𝑦, 𝑧) ∈ ℝ 3 are sufficiently smooth functions and 𝜙 and 𝜖 are rate constants that govern the separation of time scales. The variable 𝑥(𝑡) represents the membrane potential and the other two variables 𝑦(𝑡) and 𝑧(𝑡) stand for the gating dynamics of the (K + ) channels and the dynamics of cytosolic Ca 2+ , respectively. We require the right-hand sides to satisfy the following conditions: /𝑏} that guarantees a region of bistability for a range of values of 𝑧 ≥ 0 (Fig. 2). C6 The parameters 𝑎 1 , 𝑏 1 ≤ 0 and 𝑘 ≥ 0 in the (one-dimensional) slow subsystem Eqn. (3) of the polynomial model are chosen W.L.O.G. such that the 𝑧-nullcline {𝑧 = (𝑠 𝑎 1 𝑥 + 𝑏 1 ) /𝑘}, intersects the Z-shaped equilibrium curve {𝑦 = 𝑥 2 and 𝑧 = ( 𝑠 𝑎 𝑥 3 -(𝑠 + 1) 𝑥 2 ) /𝑏} of the fast subsystem Eqns. ( 1)-(2) somewhere in the middle branch, which is of saddle type (Fig. 2). This intersection point corresponds to a (degenerate) branching bifurcation of the full system Eqns. (1)-(3) at 𝜖 = 0; the bifurcation is degenerate because the Z-shaped equilibrium curve consists entirely of equilibria that have a zero eigenvalue. The branching bifurcation determines the location and stability of the fixed point FP = (𝑥 FP , 𝑦 FP , 𝑧 FP ) of the full system that persists for 𝜖 > 0; see Appendix A.2 for details.

C1 The function 𝑓 (𝑥, 𝑦, 𝑧) = -𝑠 (-𝑎 𝑥 3 + 𝑥 2 ) -𝑦 -𝑏 𝑧
C7 Plateau bursting also relies on the existence of a Hopf bifurcation in the fast subsystem (Eqns. (1)-( 2)); we assume that this Hopf bifurcation is unique.

A c c e p t e d m a n u s c r i p t

The parameter 𝑠 < 0 in the polynomial model Eqns. ( 1)-(3) plays the same role as 𝑉 m L in the generic endocrine model Eqns. ( 4)- [START_REF] Belan | Calcium clamp in isolated neurones of the snail Helix pomatia[END_REF]; it controls the type of bursting by converting the Hopf bifurcation (HB) of the fast subsystem from supercritical to subcritical (Fig. 2). The type of HB is determined by the sign of a first Lyapunov coefficient evaluated at the critical equilibrium (𝑥 HB , 𝑦 HB ) = (𝑥 HB , 𝑥 2 HB ), and in our case [START_REF] Izhikevich | Dynamical systems in neuroscience : the geometry of excitability and bursting[END_REF] it is given by:

sign [𝑙 1 (𝑥 HB )] = sign [ 𝐹 ′′′ + 𝐹 ′′ (𝐹 ′′ -𝐺 ′′ ) (𝐺 ′ -𝜙) ] = sign [ 6 𝑠 𝑎 + (2 𝑠 (3 𝑎 𝑥 HB -1)) 2 -2 (2 𝑠 (3 𝑎 𝑥 HB -1)) 2 𝑥 HB -1 ] ,
where 𝐹 (𝑥) = -𝑠 (-𝑎 𝑥 3 + 𝑥 2 ) -𝑏 𝑧 and 𝐺(𝑥) = 𝑥 2 . The values of 𝑠 < 0 that we consider are chosen W.L.O.G. such that this transition occurs in the region of bistability with respect to 𝑧 (Fig. 2).

In the following analysis we fix all the parameters in the model except for 𝜖, 𝑠 and 𝑏 1 . The parameter 𝜖 controls the speed of the slow variable 𝑧 and is our main bifurcation parameter corresponding to 𝑓 c in the generic endocrine model. The parameter 𝑠 controls the location and type of HB in the fast subsystem, which is also related to the position of the HC, in analogy with the effect that decreasing 𝑉 m L has on the behavior of the generic endocrine system (Fig. 1). Note that, similar to 𝑉 m L , the parameter 𝑠 also appears in the slow (𝑧) equation of the polynomial model. The parameter 𝑏 1 determines the location of the fixed point FP of the full system, which exists for all 𝜖 and is given by the intersection of the 𝑧-nullcline and the Z-shaped equilibrium curve of the fast subsystem Eqns. ( 1)-( 2); the locus of FP affects the bifurcations of the full system Eqns. ( 1)-(3) that occur when 𝜖 is varied. Without loss of generality and according to conditions C1-C7 we choose the rest of the system parameters to be: 𝑎 = 0.5, 𝑏 = 1, 𝑎 1 = -0.1 and 𝑘 = 0.2.

We plot in Fig. 2 the bifurcation diagram of the polynomial fast subsystem Eqns. ( 1)-( 2) using the slow variable 𝑧 as bifurcation parameter; panel (a) shows the bifurcation diagram for 𝑠 = -1.61 and panel (b) for 𝑠 = -2.6, which correspond to square-wave and pseudo-plateau bursting, respectively. A comparison between Figs. 1 and2 demonstrates that the polynomial model reproduces qualitatively the dynamics of the generic endocrine model. Similar to the biophysical system, the transition from supercritical to subcritical Hopf bifurcations in the fast subsystem of the phenomenological model Eqns. ( 1)-( 3) is accompanied by a right shift of the Z-shaped equilibrium curve that, consequently, covers a larger range of 𝑧-values during plateau bursting in the full system. In both panels of Fig. 2 we also plot three representative choices of the 𝑧-nullclines where 𝑏 1 is chosen such that these straight lines intersect the Z-curve below (green), near (blue), and well above (purple) the homoclinic bifurcation (HC) of the fast subsystem.

A c c e p t e d m a n u s c r i p t Bifurcation Analysis of the Full System

The equilibria and periodic orbits found in the polynomial fast subsystem Eqns. ( 1)-( 2) with 𝜖 = 0 provide predictions of the dynamics for 0 < 𝜖 ≪ 1 [START_REF] Bertram | Topological and phenomenological classification of bursting oscillations[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Golubitsky | An unfolding theory approach to bursting in fast-slow systems[END_REF][START_REF] Izhikevich | Neural excitability, spiking and bursting[END_REF][START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF]. We perform a numerical continuation study of the full system Eqns. ( 1)-( 3) and investigate how the periodic orbits of the full system are organized for a much larger range of 𝜖 > 0. We consider both square-wave and pseudo-plateau bursting, for 𝑠 = -1.61 and 𝑠 = -2.6, respectively (Figs. 3456).

We find that only one of the equilibria on the Z-shaped equilibrium curve persists for 𝜖 > 0, namely the fixed point FP at which the 𝑧-nullcline intersects the bifurcation diagram of the fast subsystem. If this intersection lies on the lower branch of stable equilibria or on the upper branch such that the corresponding equilibrium in the fast subsystem is stable then the full system does not support any bursting or spiking solutions. Indeed, for these cases FP is a globally stable fixed point for 𝜖 > 0 (Appendix A.2). However, if FP for 𝜖 = 0 corresponds to an unstable equilibrium of the fast subsystem, then a small perturbation 𝜖 > 0 may give rise to a periodic orbit of the full system, which corresponds to a bursting or spiking orbit [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF]. The existence and location of FP do not depend on 𝜖. However, FP eventually becomes stable for 𝜖 large enough, which is marked by a Hopf bifurcation (HB 2 ). The emanating branch of periodic orbits of the full system gives rise to a sequence of spike-adding bifurcations. The nature of this sequence appears to be determined only by the location of FP relative to the homoclinic bifurcation HC of the fast subsystem. We illustrate this in the following sections by computing the bifurcation diagram of the full model Eqns. ( 1)-( 3) for decreasing values of 𝑏 1 , which shift the locus of FP below and above the HC (Fig. 2). We then compare our findings with theoretical predictions of the behavior in the limit 𝜖 → 0. We end with a twoparameter study, where we continue the Hopf bifurcation HB 2 of the full system in the (𝜖, 𝑏 1 )-plane.

Route to Bursting via Spike-Adding Saddle-Node of Periodics Bifurcations.

We start our analysis of the bifurcation structure of the polynomial system Eqns. ( 1)-( 3) by considering the case 𝑏 1 = -0.01, for which FP lies below HC. Figure 3 shows the bifurcation diagram of the full system for 𝑠 = -1.61 and Fig. 4 for 𝑠 = -2.6. The bifurcation diagrams are presented in three-dimensional (𝜖, 𝑧, 𝑥)space with 𝜖 as the bifurcation parameter plotted on a logarithmic scale in panels (a) and on an (enlarged) linear scale in panels (b). We also plot the bifurcation diagrams of the fast subsystems Eqns. ( 1)-( 2) for 𝜖 = 0. Since on the logarithmic scale these would be pushed off to -∞, we project them onto the (𝑧, 𝑥)-plane at an arbitrary fixed value of 𝜖 = 10 -4 .

A c c e p t e d m a n u s c r i p t

The unique fixed point FP of the full system Eqns. ( 1)-( 3), which has real eigenvalues 𝜆 1 < 0 < 𝜆 2 and 𝜆 3 = 0 at 𝜖 = 0, becomes a hyperbolic saddle with 𝜆 3 > 0 for 0 < 𝜖 ≪ 1 (Appendix A.2). As 𝜖 increases, FP becomes stable in a Hopf bifurcation (HB 2 ). The Hopf bifurcation HB 2 is subcritical for both 𝑠 = -1.61 (Fig. 3) and 𝑠 = -2.6 (Fig. 4) and gives rise to a branch of unstable periodic orbits that becomes stable in a saddle-node of periodics (SNP 1 ). The branch of stable periodic orbits corresponds to tonic spiking of large amplitude, unlike the tonic spiking typically seen in square-wave bursting models as the Ca 2+ pump rate is increased [START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF] or the conductance of K Ca channels is decreased [START_REF] Chay | On the effect of the intracellular calcium-sensitive K + channel in the bursting pancreatic 𝛽-cell[END_REF][START_REF] Medvedev | Reduction of a model of an excitable cell to a onedimensional map[END_REF]. This branch can be considered as bursts with one spike. Sample spiking trajectories of the full system are superimposed on the bifurcation diagrams in Figs. 3(b) and 4(b). The first one (from the right) in both figures is a two-spike periodic orbit. As 𝜖 decreases, the branch of one-spike periodic orbits loses and regains stability in saddle-node of periodics (SNP) bifurcations and during this process it transforms from a one-spike into a two-spike periodic orbit. The transition happens over a very narrow range of 𝜖, to the right of SNP 2 , during which the stable one-spike periodic orbit coexists with a stable two-spike periodic orbit. As 𝜖 decreases further the series of SNP bifurcations repeats, delimiting smaller and smaller portions of the branch, each of which corresponds to a bursting solution with one more spike (Figs. 3 and4). Using the software package AUTO [START_REF] Doedel | AUTO97: continuation and bifurcation software for ordinary differential equations[END_REF], we were able to follow this branch down to 𝜖 = 𝑂(10 -3 ). Figures 3(b) and 4(b) illustrate the accumulation of the SNP bifurcations as 𝜖 decreases for square-wave (𝑠 = -1.61) and pseudoplateau (𝑠 = -2.6) bursting, respectively. The bifurcation diagrams for squarewave and pseudo-plateau bursting are very similar and both exhibit a sequence of SNP bifurcations creating 𝑛-spike solutions for increasingly larger 𝑛 as 𝜖 decreases. In the limit 𝜖 → 0 the number of spikes of the stable bursting solutions goes to infinity, while the stability region of each individual orbit goes to zero.

Figures 3 and4 suggest that the full system exhibits a spike-adding cascade (2, 7, 35, 52) mediated by SNP bifurcations if the fixed point FP of the full system lies below the homoclinic bifurcation HC of the fast subsystem for 𝜖 = 0 (Figs. 2(a) and (b) with 𝑏 1 = -0.01).

Route to Bursting via Spike-Adding Isolas.

We continue our analysis of the bifurcation structure of the polynomial model Eqns. ( 1)-( 3) by considering the cases 𝑏 1 = -0.045, 𝑠 = -1.61 and 𝑏 1 = -0.21, 𝑠 = -2.6, for which FP lies above HC (Figs. 5 and6, respectively). As before, we use 𝜖 as the bifurcation parameter and plot the bifurcation diagrams of the full system in (𝜖, 𝑧, 𝑥)-space on a logarithmic scale in panels (a) and on an (enlarged) linear scale in panels (b). The bifurcation diagrams of the fast subsystems Eqns. Our choices for 𝑏 1 and 𝑠 illustrate the two possibilities for positioning the unique fixed point FP of the full system Eqns. ( 1)-(3). For 𝑏 1 = -0.045 and 𝑠 = -1.61 the 𝑧-nullcline intersects the bifurcation diagram of the fast subsystem above HC, but below SN 1 . For 𝑏 1 = -0.21 and 𝑠 = -2.6, on the other hand, this intersection lies between SN 1 and HB 1 ; compare Figs. 2(a) and (b). As before, if FP lies between HC and SN 1 , it has real eigenvalues 𝜆 1 < 0 < 𝜆 2 and 𝜆 3 = 0 at 𝜖 = 0 and becomes a hyperbolic saddle with 𝜆 3 > 0 for 0 < 𝜖 ≪ 1 (Appendix A.2). However, if FP lies between SN 1 and HB 1 , the eigenvalues from the fast subsystem are unstable; they are real for FP close to SN 1 and form a complex conjugate pair for FP close to HB 1 . In this case 𝜆 3 < 0 and FP is again a hyperbolic saddle with two unstable eigenvalues for 0 < 𝜖 ≪ 1 (Appendix A.2).

We find that the bifurcation diagram of the full system is topologically equivalent for these two choices of FP above HC, but rather different from the case where FP lies below HC. Note that, locally near 𝜖 = 0, there is no difference in whether FP lies below or above HC; as 𝜖 increases FP becomes a hyperbolic saddle with two unstable eigenvalues and it gains stability in a Hopf bifurcation (HB 2 ). However, in contrast to Figs. 3 and4, the Hopf bifurcation HB 2 is now supercritical and gives rise to a branch of stable periodic orbits that correspond to large-amplitude tonic spiking. As before, we superimpose sample spiking trajectories of the full system on the bifurcation diagrams in Figs. 5(b) and 6(b) and an example of a stable one-spike periodic orbit is shown in Fig. 6(b), where 𝑠 = -2.6.

For both 𝑠 = -1.61 and 𝑠 = -2.6 the branch of one-spike periodic orbits loses stability in a supercritical period-doubling bifurcation (PD 1 ). The emanating branch of stable period-doubled orbits corresponds to two-spike periodic orbits, examples of which are superimposed on the bifurcation diagrams in Figs. 5(b) and 6(b). The two-spike periodic orbit loses stability in another period-doubling bifurcation (PD 2 ) that gives rise to a period-doubled two-spike orbit; an example of such a periodic orbit is shown in Fig. 5(b) and it does not correspond to a standard bursting solution. The period-doubled two-spike orbit is stable for a much smaller range in 𝜖 and it also loses stability in a period-doubling bifurcation, starting what appears to be a period-doubling cascade (not shown). We refer to these and further period-doubled 𝑛-spike orbits as secondary bursting solutions and we do not pursue further investigation of these types of period-doubled orbits.

Instead, we focus on the spike-adding cascade that occurs also if FP lies above HC. In this case, the spike-adding cascade is organized by isolas. A family of stable 𝑛-spike periodic orbits is born in an SNP bifurcation and, as 𝜖 decreases, it undergoes period-doubling bifurcations to secondary bursting solutions. We found these isolas by generating seed solutions using numerical integration of Eqns. (1)-

A c c e p t e d m a n u s c r i p t

(3) for decreasing fixed values of 𝜖 and continuing them in AUTO [START_REF] Doedel | AUTO97: continuation and bifurcation software for ordinary differential equations[END_REF]. Figures 5 and6 show isolas of 𝑛-spike periodic orbits for 𝑛 = 3, . . . , 9 and 𝑛 = 3, . . . , 12, respectively, along with examples of bursting trajectories with increasing numbers of spikes. We observe that the isolas create gaps between stable 𝑛-spike periodic orbits for moderately small values of 𝜖, but they overlap as 𝜖 decreases creating narrow intervals with coexisting stable 𝑛and (𝑛 + 1)-spike periodic orbits. However, the stable portions of the isolas become smaller as 𝑛 increases resulting in smaller regions in parameter space where each stable 𝑛-spike periodic orbits exist. Figures 5 and6 suggest that the full system exhibits a spike-adding cascade mediated by period-doubling bifurcations and isolas if the fixed point FP of the full system lies above the homoclinic bifurcation HC of the fast subsystem for 𝜖 = 0 (Figs. 2(a) and (b) with 𝑏 1 = -0.045 and 𝑏 1 = -0.21, respectively). In particular, it does not seem to matter whether FP lies below or above the saddlenode bifurcation SN 1 of the fast subsystem for 𝜖 = 0 as long as the corresponding equilibrium of the fast subsystem is unstable.

Chaotic Bursting Solutions

When the route to bursting is mediated via spikeadding SNP bifurcations (FP well below HC), periodic bursting is accompanied by bistability and chaotic alternation between regular 𝑛and (𝑛 + 1)-spike periodic orbits. This has previously been shown for fold-homoclinic bursting [START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF], and an example of irregular, presumably chaotic, alternation between two-and three-spike solutions is illustrated for fold-subHopf bursting in Fig. 7(a). In contrast, when bursting arises via spike-adding isolas (FP well above HC), bursting can be chaotic due to the overlapping of isolas in regimes where period-doubling cascades exist. These give rise to chaotic alternations between regular 𝑛-spike periodic orbits and secondary bursting solutions. As an illustration of such behavior, we plot in Fig. 7(b) a time series of the polynomial model for the fold-subHopf case, showing a spontaneous transition from period-doubled two-spike solutions to three-spike bursting.

Mixed Route to Bursting. When FP lies close to the HC for 𝜖 = 0, the periodic solutions branches are of mixed type that is characterized by the presence of both PDs and SNPs in the bifurcation diagram of the full system. As an illustration we computed the bifurcation diagrams of the full system Eqns. ( 1)-( 3) both for the cases of fold-homoclinic bursting, with 𝑠 = -1.61 and 𝑏 1 = -0.023 (Fig. 2(a)), and of fold-subHopf bursting, with 𝑠 = -2.6 and 𝑏 1 = -0.066 (Fig. 2(b)).

The fold-homoclinic case is shown in Fig. 8(a) and corresponds to a situation where FP lies just below the HC for 𝜖 = 0. In contrast to Fig. 3, the Hopf bifurcation HB 2 is supercritical and gives rise to a branch of stable one-spike periodic

A c c e p t e d m a n u s c r i p t

orbits that ends at a period-doubling bifurcation (PD 1 ). For these relatively large values of 𝜖, the bifurcation diagram resembles that of Fig. 5: a branch of stable two-spike periodic orbits emanates from PD 1 that loses stability in another perioddoubling bifurcation (PD 2 ), which gives rise to secondary bursting solutions. In a narrow 𝜖-interval these secondary bursting solutions coexist with a branch of stable three-spike periodic orbits. As for Fig. 5, we did not further investigate the period-doubling cascade of secondary bursting solutions, but rather concentrated on the spike adding. The branch of three-spike periodic orbits (magenta curve) in Fig. 8(a) does not lie on an isola, in contrast to Fig. 5. For these smaller values of 𝜖 the bifurcation diagram resembles that of Fig. 3, as expected. Instead of individual isolas, continuation of the three-spike periodic orbits leads to a single connected branch of 𝑛-spike periodic orbits that consists of increasing numbers of spikes as 𝜖 decreases. We were able to follow this branch down to values of 𝜖 = 𝑂(10 -3 ).

The fold-subHopf case is shown in Fig. 8(b) and corresponds to a situation where FP lies slightly above HC for 𝜖 = 0. Here, HB 2 is also supercritical and the stable branch of one-spike periodic orbits again loses stability in a period-doubling bifurcation (PD). However, the emanating branch of stable two-spike periodic orbits undergoes a sequence of SNP bifurcations corresponding to a spike-adding cascade and the entire family of periodic orbits in the full system forms a single connected branch. As before we were only able to follow the branch down to values of 𝜖 = 𝑂(10 -3 ).

The above computations indicate that there is an interesting transition between the two routes to bursting in both classes of models as FP crosses from one side of HC to the other. Detailed investigation of this transition is left for future investigation, because it requires numerical exploration in a region of very small values of 𝜖 where our computations break down.

Behaviour in the limit of small 𝜖.

We have studied the case of 𝜖 = 0 (fast subsystem bifurcations) and the cascade of periodic orbit bifurcations as 𝜖 decreases from large values, but it is evident that there are important phenomena in the region of small 𝜖 that our numerical continuations have not addressed. We extend here the theory in [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF] for 0 < 𝜖 ≪ 1, which only applies to the case where FP is located just below or above the HC at 𝜖 = 0, and discuss how it ties in with our numerical study for a much larger range of 𝜖.

The limit of small 𝜖 for square-wave bursting Let us first consider square-wave (fold-homoclinic) bursting, which is Scenario 3 in [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF]. Terman [START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF] considered this case already in 1992, but we will follow Belykh et al. [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF]. If FP lies just below the

A c c e p t e d m a n u s c r i p t

HC (the case 𝑙 < 0 in (7)), then there exists a bursting solution for 0 < 𝜖 ≪ 1 small enough. In the limit 𝜖 → 0 this bursting solution accumulates on a periodic orbit with infinitely many spikes and the range of the slow variable (𝑧 in our case) covers the interval between the lower saddle-node bifurcation SN 2 and the HC [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF].

Our numerical study indicates that the theoretical predictions for 0 < 𝜖 ≪ 1 are also valid when FP lies much closer to SN 2 . The continuations shown in Figs. 3 and8(a) indicate that the predicted branch of bursting solutions connects via a single cascading family of 𝑛-spike periodic orbits to the family of periodic orbits that emanates from the Hopf bifurcation HB 2 at large 𝜖.

If FP lies just above the HC then Belykh et al. [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF] predict the existence of bursting solutions, but the singular limit corresponds to a regime where continuous (tonic) spiking exists. This tonic spiking is a periodic orbit with an amplitude that is close to the homoclinic orbit at HC for 𝜖 = 0 and is different from the branch of large-amplitude tonic spiking solutions that we found emanating from the Hopf bifurcation HB 2 at large 𝜖. The periodic orbit for 0 < 𝜖 ≪ 1 may lose stability in a period-doubling bifurcation, and it definitely does not persist beyond an SNP bifurcation that is predicted to occur for some value 𝜖 ≪ 1. We were unable to identify numerically the end of the regime of tonic spiking. The bursting oscillations arise from a homoclinic bifurcation of the full system, where the one-dimensional stable manifold of FP is contained in its two-dimensional unstable manifold. This homoclinic bifurcation happens before the SNP bifurcation predicted by Belykh et al. [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF], at a value of 0 < 𝜖 ≪ 1 that is expected to depend on the distance between FP and the HC at 𝜖 = 0. Hence, the regimes of tonic spiking and bursting may overlap, but they do not form a single connected family. We have not been able to locate the homoclinic bifurcation of FP in the full system, but our numerical findings indicate that the bursting regime is organized by a family of spike-adding isolas that starts with the family of one-spike periodic orbits emanating from HB 2 at relatively large 𝜖 and accumulates on this homoclinic bifurcation of FP at 0 < 𝜖 ≪ 1 (Fig. 5).

It is again possible to extend the theory for 0 < 𝜖 ≪ 1 to cases where FP lies well above HC. The periodic orbit that corresponds to tonic spiking emanates from a closed curve at 𝜖 = 0 that is some sort of average of a small family of the stable periodic orbits in the fast subsystem, where the range of the family depends on where the 𝑧-nullcline intersects this family; see [START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF]. Note that this argument also holds in the range for 𝑏 1 such that FP lies between HB 1 and SN 1 at 𝜖 = 0. As FP is moved further away above the HC, we expect that continuous spiking persists for increasingly larger values of 0 < 𝜖 ≪ 1. Indeed, both Terman (46) and Belykh et al. [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF] predict that, for fixed 0 < 𝜖 ≪ 1, a continuous variation from FP above the HC to FP below the HC (in our case this means increasing the parameter 𝑏 1 ) leads to a transition from continuous spiking to bursting via a regime with chaotic dynamics; an unpredictable number of spikes within each burst occurs in regimes

A c c e p t e d m a n u s c r i p t

where the isolas overlap (52).

The limit of small 𝜖 for pseudo-plateau bursting

There is yet no predictive theory for the range 0 < 𝜖 ≪ 1 for pseudo-plateau (fold-subHopf) bursting. Pseudo-plateau bursting in the biophysical sense (small transient spikes) relies on the fact that 𝜖 is only moderately small, but the fold-subHopf structure persists for 0 < 𝜖 ≪ 1. For such very small values of 𝜖, the case of fold-subHopf bursting compares to Scenario 1 in [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF], that is, for 0 < 𝜖 ≪ 1 the solutions are relaxation oscillations determined by the branches of stable equilibria; for Scenario 1 the two branches of stable equilibria end at the saddle-node bifurcations SN 1 and SN 2 , but in our case the upper equilibrium branch loses stability already at the subcritical Hopf bifurcation HB 1 . Hence, we expect the scenario of slow passage through a Hopf bifurcation (3) and the bursting solutions for 0 < 𝜖 ≪ 1 to resemble relaxation oscillations that exhibit many oscillations during the depolarized-plateau phase with amplitudes that first decrease exponentially and subsequently increase exponentially.

We can get a glimpse into the region with 0 < 𝜖 ≪ 1 via selected numerical integrations in time. Inspection of the slow-variable (𝑧) oscillations for fold-subHopf bursting suggests persistence of slow oscillations with periods going to infinity. The slow passage through the Hopf bifurcation HB 1 anticipated above, results in the maximum of 𝑧 occurring at the same distance from the 𝑧-value of HB 1 as its minimum. If HB 1 lies too close to SN 1 then the maximum of 𝑧 equals the 𝑧-value of SN 1 . Figure 9(a) shows a bursting solution for the polynomial model with 𝑠 = -2.6, 𝑏 1 = -0.01 and 𝜖 = 10 -4 . The bursting solution resembles a relaxation oscilation, but the depolarized-plateau phase is characterized by a slow passage through the Hopf bifurcation; note that the pseudo spikes of the 𝑥-variable during the depolarized-plateau phase have nearly disappeared, with only vestigial spikes at the beginning and end of the plateau. For some values of 𝜖 and 𝑏 1 , the pseudo-plateau bursting appears to be chaotic, exhibiting depolarized plateaus with variable and unpredictable duration. That is, the minimum of 𝑧 lies at SN 2 , but its maxima lie between HB 1 and SN 1 . An example is shown in Fig. 9(b) with 𝑠 = -2.6, 𝑏 1 = -0.12 and 𝜖 = 0.001.

The Hopf bifurcation of the full system

In all the examples that we showed of the bifurcation diagam of the full system Eqns. ( 1)-( 3), the Hopf bifurcation HB 2 of the fixed point FP happens at a relatively large value of 𝜖. However, depending on the choice for 𝑏 1 , which moves the location of FP on the Z-shaped equilibrium curve at 𝜖 = 0 relative to the HC of the fast subsystem, HB 2 can occur for arbitrarily small 0 < 𝜖 ≪ 1. If we increase 𝑏 1 to 0 (this means going down in Fig. 10), starting from HB 2 at 𝑏 1 = -0.01 where FP lies below the HC at 𝜖 = 0 for both choices of 𝑠 (Figs. 3 and4), the Hopf bifurcation curve ends at the saddle-node bifurcation SN 2 at 𝜖 = 0. As discussed in Appendix A.2, the point SN 2 at 𝜖 = 0 corresponds to a singular Hopf bifurcation of the full system, that is, the two-dimensional center manifold of the Hopf bifurcation involves both slow and fast directions. The singular Hopf bifurcation persists for 0 < 𝜖 ≪ 1, where it occurs at a value of 𝑏 1 for which FP lies 𝑂(𝜖) close to (but past) the left knee of the Z-shaped equilibrium curve. The occurrence of a singular Hopf bifurcation has been shown previously for the original Hindmarsh-Rose system (26) in (2, Fig. 12) as well as in [START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF]; see also [START_REF] Guckenheimer | Singular Hopf bifurcation in systems with two slow variables[END_REF]. Note that HB 2 is only a singular Hopf bifurcation for 0 < 𝜖 ≪ 1 and -1 ≪ 𝑏 1 ≤ 0, but the exact transition from an ordinary to a singular Hopf bifurcation is not defined. When FP corresponds to a saddle very close to SN 2 at 𝜖 = 0, it is unstable for 0 < 𝜖 ≪ 1, but its two unstable eigenvalues are complex conjugate and lie extremely close to the imaginary axis. Hence, only a very small increase in 𝜖 already stabilizes FP as the singular Hopf bifurcation HB 2 occurs. The singular Hopf bifurcation gives rise to a small-amplitude periodic orbit that transforms very quickly as 𝜖 varies over an exponentially small interval; if HB 2 is supercritical, 𝜖 will be decreasing, but if HB 2 is subcritical the branch will turn around at an SNP bifurcation that happens exponentially close after HB 2 . For the fold-homoclinic case, the periodic orbit transforms into an 𝑛-spike bursting orbit where 𝑛 is extremely large for 𝜖 small. For the fold-subHopf case, the periodic oribt transforms into a relaxation oscillation. These exponentially small transitions involve so-called canard orbits, where the periodic orbits contain segments that follow the saddle-branch of the Z-shaped steady-state curve [START_REF] Guckenheimer | Singular Hopf bifurcation in systems with two slow variables[END_REF].

If we decrease 𝑏 1 from 𝑏 1 = -0.01, we find that the Hopf curve for both 𝑠 = -1.61 and 𝑠 = -2.6 initially increases in 𝜖, but then decreases again until it ends at 𝜖 = 0. However, the Hopf curve is not monotonic in 𝑏 1 . Initially, the fixed point FP moves up the middle branch at 𝜖 = 0 and past SN 1 onto the upper branch of unstable equilibria until it reaches HB 1 at a value of 𝜖 > 0. However, the Hopf curve continues for decreasing 𝜖 when we trace FP past HB 1 onto the stable segment of the upper branch of the Z-shaped curve. As 𝜖 decreases further, the Hopf curve reaches a minimum in 𝑏 1 (which corresponds to a maximum in Fig. 10) and then returns to the value of 𝑏 1 that corresponds to HB 1 as 𝜖 → 0. Hence, for a small range of 𝑏 1 when FP is located close to but to the left of HB 1 on the stable segment of the upper equilibrium branch, the point FP is initially stable and there exists a small range of values 0 < 𝜖 ≪ 1 for which FP is unstable and the attractor is a small-amplitude periodic orbit; the end points of this small interval are

A c c e p t e d m a n u s c r i p t

(ordinary) Hopf bifurcation points. When FP is located on the unstable segment of the upper branch, then it is initially unstable and only one Hopf bifurcation exists for 0 < 𝜖 ≪ 1, which stabilizes FP as 𝜖 increases.

Note that in the fold-homoclinic case HB 1 is supercritical, while in the fold-subHopf case it is subcritical. Furthermore, HB 2 is subcritical when FP at 𝜖 = 0 lies near SN 2 (see Fig. 3), but it is supercritical when FP at 𝜖 = 0 lies near SN 1 (see Fig. 6). Hence, there must be at least one degenerate Hopf bifurcation on the blue curve (𝑠 = -1.61) corresponding to the fold-homoclinic case, and at least two degenerate Hopf bifurcations on the red curve (𝑠 = -2.6) corresponding to the fold-subHopf case in Fig. 10.

Discussion

Given the importance of bursting activity in endocrine cells, which sets the levels of [Ca 2+ ] i [START_REF] Beauvois | Glucose-induced mixed [Ca 2+ ] c oscillations in mouse betacells are controlled by the membrane potential and the SERCA3 Ca 2+ -ATPase of the endoplasmic reticulum[END_REF][START_REF] Bergsten | Role of Oscillations in Membrane Potential, Cytoplasmic Ca 2+ , and Metabolism for Plasma Insulin Oscillations[END_REF][START_REF] Van Goor | Dependence of pituitary hormone secretion on the pattern of spontaneous voltage-gated calcium influx. cell type-specific action potential secretion coupling[END_REF][START_REF] Vincent | Endocrine cell excitability opens the way to novel pharmacological intervention: example of the anterior pituitary cell[END_REF] and is instrumental for the regulation of hormone exocytosis, it is of interest to identify the key mechanisms governing it. We considered two general classes of models for endocrine bursting, square-wave bursting models (fold-homoclinic) and pseudo-plateau bursting models (fold-subHopf). It is important to understand the similarities and differences between square-wave and pseudo-plateau bursting because different kinds of endocrine cells may be tuned to operate in one or the other regime depending on their functions. For example, the models indicate that pseudo-plateau bursting may provide higher Ca 2+ for secretion. We have presented here a generic Hodgkin-Huxley type model that captures the main features of a number of previously published models and used it to identify variation in properties of the inward, excitatory current, specifically the parameter 𝑉 m L , as a possible way to switch between the classes (see below). We also used a simplified polynomial model [START_REF] Hindmarsh | A model of neuronal bursting using 3 coupled 1st order differential-equations[END_REF], which emphasized the general dynamic features of the two classes of bursters and was more convenient for the challenging numerical continuations carried out here.

Physiological implications.

Although we have focused on the mathematical effects of varying the parameter 𝑓 c (corresponding to 𝜖 in the simplified polynomial model), it represents the fraction of free cytosolic Ca 2+ and accounts for the buffering capacity of cells. Smaller values of 𝑓 c slow down the rise in Ca 2+ and in turn the activation of the K Ca channels. As shown here, when 𝑓 c is small the models exhibit bursts with more spikes and, hence, longer depolarized plateaus. This agrees with a recent study [START_REF] Roussel | Modulation of neuronal excitability by intracellular calcium buffering: from spiking to bursting[END_REF] that combined modeling and experiments to show that cytosolic calcium buffering capacity can tightly modulate neuronal firing patterns and determine whether bursting or spiking is generated. The range

A c c e p t e d m a n u s c r i p t

of parameter values (𝑓 c ≥ 10 -3 , 𝜖 ≥ 10 -3 ) we were able to explore numerically is comparable to that found in most previously published models (ranging from 10 -3 to 10 -1 ) [START_REF] Bertram | The phantom burster model for pancreatic beta-cells[END_REF][START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF][START_REF] Chay | Minimal model for membrane oscillations in the pancreatic beta-cell[END_REF][START_REF] Lebeau | Analysis of a reduced model of cortiocotroph action potentials[END_REF][START_REF] Stern | Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus[END_REF][START_REF] Tabak | Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents[END_REF][START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF][START_REF] Tsaneva-Atanasova | Diffusion of calcium and metabolites in pancreatic islets: killing oscillations with a pitchfork[END_REF][START_REF] Van Goor | Paradoxical role of large-conductance calcium-activated K + (BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells[END_REF] and also observed in cells [START_REF] Al-Baldawi | Cytoplasmic calcium buffer capacity determined with Nitr-5 and DM-nitrophen[END_REF][START_REF] Belan | Calcium clamp in isolated neurones of the snail Helix pomatia[END_REF][START_REF] Berlin | Intrinsic cytosolic calcium buffering properties of single rat cardiac myocytes[END_REF][START_REF] Fierro | High endogenous calcium buffering in Purkinje cells from rat cerebellar slices[END_REF][START_REF] Fleet | Calcium buffering capacity of neuronal cell cytosol measured by flash photolysis of calcium buffer NP-EGTA[END_REF][START_REF] Müller | Calcium buffering in bursting Helix pacemaker neurons[END_REF][START_REF] Neher | Calcium gradients and buffers in bovine chromaffin cells[END_REF][START_REF] Schwiening | Relationship between intracellular calcium and its muffling measured by calcium iontophoresis in snail neurones[END_REF][START_REF] Zhou | Mobile and immobile calcium buffers in bovine adrenal chromaffin cells[END_REF]. Calcium buffering capacity is not only variable among cells but can change under different physiological conditions, such as the developmental stage. In a recent study in hippocampal granule cells, younger cells had approximately three times smaller Ca 2+ -binding ratio (43) than older cells. We note that the unusually slow Ca 2+ oscillations in pancreatic 𝛽-cells (periods from tens of seconds to several minutes) are likely not the result of very small binding ratio but rather reflect the slow dynamics of metabolism and/or of Ca 2+ in the endoplasmic reticulum [START_REF] Bertram | A calcium-based phantom bursting model for pancreatic islets[END_REF].

A key difference between the two types of bursters is that the spikes disappear in the fold-subHopf case as 𝑓 c or 𝜖 goes to 0. It may thus be possible to distinguish the two types experimentally by reducing 𝜖 via addition of exogenous Ca 2+ buffer as in [START_REF] Roussel | Modulation of neuronal excitability by intracellular calcium buffering: from spiking to bursting[END_REF]. This also sheds new light on a question raised by the early observation that square-wave bursters are like relaxation oscillators with spikes superimposed on the upper state. Namely, if cells need a plateau to raise [Ca 2+ ] i , why do they not just have a plateau instead of spikes? One possibility is that cells have to build the plateau out of ion channel interactions, and the delayed rectifier is too slow to cancel out the spikes and produce a pure plateau. In addition, given the naturally occurring range of 𝑓 c (10 -3 -10 -1 ), Ca 2+ is too fast compared to the rate of attraction to the upper state to eliminate the spikes. With the earlier foldhomoclinic models, this was not so apparent because the spike amplitude does not vary with 𝑓 c . In this sense, one can view fold-subHopf as a case that is closer to a plateau and in fact results in higher calcium, at least when mediated by a shift in 𝑉 m L . It is also possible that the deeper spike repolarizations found in square-wave bursting slow down or reverse the rise of Ca 2+ and thus prolong the active phases of the bursts.

We have identified one route to transform fold-homoclinic to fold-subHopf bursting, namely, translation of the activation curve of the Ca 2+ current; there may be other parameters that can achieve this. Changing the conductance of voltagegated K + channels or K Ca channels, such as the BK channel, which is voltageas well as Ca 2+ -sensitive, can also shift the slow manifolds, change spike amplitude, and change the location and nature of the fast-subsystem Hopf bifurcations in endocrine cell models. BK blockade and natural variation of BK channel density were shown to affect the period of apparent pseudo-plateau bursting in pituitary somatotrophs [START_REF] Tsaneva-Atanasova | Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory[END_REF], and BK blockade was shown to convert apparent bursting to large amplitude spiking [START_REF] Van Goor | Paradoxical role of large-conductance calcium-activated K + (BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells[END_REF]. Similar effects are seen by varying the time constant of voltage-gated K + channels or the conductance of voltage-dependent Ca 2+ channels (unpublished observations). None of these changes has been seen to convert

A c c e p t e d m a n u s c r i p t

fold-subHopf to fold-homoclinic bursting or vice versa, though it is possible that they could do so in combination with other parameter changes.

Full-system bifurcation analysis and spike adding. In this paper we have examined both the square-wave and pseudo-plateau bursting regimes in terms of bifurcations of the full polynomial system over a relatively large range of 𝜖. Similar analyses have been carried out previously for fold-homoclinic bursting by fixing 𝜖 and varying other parameters (2, 7, 52), but we have treated fold-subHopf bursting for the first time. We considered the two together as they form a class of closely related biological systems. We showed the emergence of bursting from a primordial large-amplitude spiking solution that undergoes a complicated cascade of bifurcations in the full system as the parameter 𝜖 decreases (Figs. 3456). Several period-doubling bifurcations that give rise to two-, four-and higher-periodic orbits and could be the start of a classic period-doubling cascade to chaotic orbits, spawn a cascade of spike-adding bifurcations that generates a sequence of new period-2,3,4,5,...,𝑛 attractors, for some finite number 𝑛, as 𝜖 decreases. The latter correspond to bursting trajectories with the respective number of spikes. The increase in the number of spikes as 𝜖 decreases has the geometric interpretation that, as the slow-variable component (the 𝑧-coordinate of the polynomial model) of the bursting orbits of the full system slows down, the trajectory spends more and more time moving along the 𝑧 direction in phase space.

Previous studies have examined spike adding for fold-homoclinic bursting [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF][START_REF] Wang | Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle[END_REF]. Here, Figs. 4, 6 and 8(b) show that spike adding happens similarly as 𝜖 is reduced in the case of fold-subHopf bursting. We showed, in addition, for both types of bursters that the form of the transition depends on where the fullsystem fixed point FP at 𝜖 = 0 is located relative to the fast-subsystem HC. When FP at 𝜖 = 0 lies well below the HC (a case that has been considered for foldhomoclinic bursters in [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable-membranes[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF][START_REF] Wang | Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle[END_REF]), the periodic branch is connected and the addition of each spike is marked by a pair of SNP bifurcations (Fig. 4). In contrast, when FP at 𝜖 = 0 lies well above the HC, the periodic branch is disjoint, with isolas and period-doubling bifurcations for each new spike (Fig. 6).

Further, as shown in Fig. 8, where FP at 𝜖 = 0 and the HC are very close, the spike-adding transition appears to be of mixed type. Initially, as 𝜖 decreases the usual period-doubling cascade occurs that characterizes the transition mediated via isolas. However, we find no individual isolas for 3,4,5,..-spike solutions but a connected complicated periodic branch/isola that comprises these solutions and features saddle-node of periodics rather than period-doubling bifurcations. This suggests that there is a continuous transition between the two mechanisms of spike

A c c e p t e d m a n u s c r i p t

adding. Detailed analysis of this transition is left for future investigations as it requires analysis of the small 𝜖-limit where our numerical computations break down. Whereas the spike-adding cascades are similar for fold-homoclinic and fold-subHopf bursters, the two differ in their behaviors as 𝜖 → 0. As shown previously [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF], when FP at 𝜖 = 0 lies above the HC, fold-homoclinic bursting must give way to small-amplitude continous spiking as 𝜖 → 0. In contrast, fold-subHopf bursters continue to exhibit large-amplitude oscillations as long as FP at 𝜖 = 0 lies to the right of HB 1 . In the fold-subHopf case, the number of spikes increases to infinity as 𝜖 → 0 as in the fold-homoclinic case, but their amplitude decreases to 0 because the rate of increase of the slow variable (𝑧 or Ca 2+ ) becomes much smaller than the rate of attraction to the branch of stable equilibria of the fast subsystem. Thus, for small 𝜖 we find plateaus with no spikes (Fig. 9(a)).

Previous studies of spike adding have been done for the fold-homoclinic case by studying 1D Poincaré return maps or specially constructed 2D maps [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Channell | Origin of bursting through homoclinic spike adding in a neuron model[END_REF][START_REF] Medvedev | Reduction of a model of an excitable cell to a onedimensional map[END_REF][START_REF] Medvedev | Transition to bursting via deterministic chaos[END_REF][START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF][START_REF] Wang | Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle[END_REF], in contrast to our approach of continuing periodic orbits of the full ODE system. Several studies have also constructed bifurcation diagrams of the full system of ODEs [START_REF] Alexander | On the dynamics of bursting systems[END_REF][START_REF] Channell | Origin of bursting through homoclinic spike adding in a neuron model[END_REF][START_REF] Duan | Two-parameter bifurcation analysis of firing activities in the Chay neuronal model[END_REF]. Such studies have pointed to the importance of spike-adding transitions for the genesis of bursting from simple continuous spiking solutions. Most studies have used a primary bifurcation parameter that translates the slow-variable nullcline; Ca 2+ pump rate or the conductance of K Ca channels have been popular choices, motivated by the ability of these parameters to convert bursting to small-amplitude continuous spiking in the first biophysical model of fold-homoclinic bursting [START_REF] Chay | Minimal model for membrane oscillations in the pancreatic beta-cell[END_REF].

Taking these studies together with ours, we can identify three distinct transitions between bursting and spiking in the fold-homoclinic case. When 𝜖 is small, there is a transition from bursting to small-amplitude spiking as FP at 𝜖 = 0 crosses the HC from below. In the limit as 𝜖 → 0, the transition occurs precisely when FP coincides with the saddle equilibrium that corresponds to the HC [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF], but bursting can persist when 𝜖 > 0 for a range of FP loci above the HC. This is the second transition: when FP at 𝜖 = 0 lies above the HC, spiking converts to bursting as 𝜖 increases from 0, passing through a complicated chaotic cascade along the way that includes plateaus of arbitrary and fluctuating duration [START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF]. Finally, as 𝜖 becomes larger, the number of spikes per burst decreases, reaching a region of large-amplitude spiking (the third transition) before oscillations end in a Hopf bifurcation (HB 2 in Figs. 3-6 and8). Those spiking solutions can be thought of as infinite trains of bursts with one spike, in contrast to the small-amplitude spiking above, which can be thought of as single bursts with an infinite number of spikes. Of the above transitions, only the third occurs in fold-subHopf bursting due to the lack of stable limit cycles in the fast subsystem.

A c c e p t e d m a n u s c r i p t

Future directions. Our analysis is complementary to the study of Belykh et al. [START_REF] Belykh | Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models[END_REF], who did a local analysis for small 𝜖; namely, we considered the behavior at large 𝜖. Fold-homoclinic bursting was addressed in their Scenario 3 and, as suggested above, fold-subHopf bursting has some similarities to their Scenario 1, though the latter is a pure relaxation oscillator with no Hopf bifurcation in the fast subsystem (HB 1 in our nomenclature). The present study leaves a gap between the very small and moderate 𝜖 regimes, and new numerical approaches to fill this gap, perhaps with appropriately rescaled equations, are suggested as a fruitful area for future study.

We have focused on full-system bifurcations here, but we also recognize that the differences between fold-homoclinic and fold-subHopf bursters lie in the fast sub-system, namely, the character of the Hopf bifurcation HB 1 . Moreover, the differences in how bursts are terminated in the two classes of systems will require analysis at the fast-subsystem level. A promising framework in which to investigate this question is that of Golubitsky et al. [START_REF] Golubitsky | An unfolding theory approach to bursting in fast-slow systems[END_REF], where local unfoldings of singularities in the fast subsystem bifurcations were used to classify types of bursting. current to flux, where 𝐴 cell = 𝜋 × 𝑑 2 cell is the surface area of the cell; and 𝑘 PMCA is the plasma membrane Ca 2+ ATPase pump rate. Since 𝑐 represent the free Ca 2+ concentration in the cytosol the corresponding fluxes in Eq. ( 6) are multiplied by the fraction 𝑓 c of free to total cytosolic Ca 2+ . The currents included in the model equations are:

𝐼 Ca (𝑉 m ) = 𝑔 Ca 𝑚 2 ∞ (𝑉 m ) (𝑉 m -𝑉 Ca ) , ( 7
) 𝐼 K (𝑉 m , 𝑛) = 𝑔 K 𝑛 (𝑉 m -𝑉 K ) , ( 8
) 𝐼 K Ca (𝑉 m , 𝑐) = 𝑔 K Ca 𝑠 ∞ (𝑐) (𝑉 m -𝑉 K ) . (9)
The steady state activation functions are:

𝑚 ∞ (𝑉 m ) = ( 1 + exp ( 𝑉 m L -𝑉 m 𝑠 𝑚 )) -1 , ( 10
)
𝑛 ∞ (𝑉 m ) = ( 1 + exp ( 𝑉 n -𝑉 m 𝑠 𝑛 )) -1 , (11)
𝑠 ∞ (𝑐) = 𝑐 4 𝑐 4 + 𝑘 4 s . (12)
Values of all parameters used in the model simulations are given in Table A.1. 1)-( 2) at 𝜖 = 0 and can be controlled by a single parameter. A convenient choice is the parameter 𝑏 1 as shown in Fig. 2. Hence, we assume that 𝑎 = 0.5, 𝑏 = 1, 𝑎 1 = -0.1, and 𝑘 = 0.2 are fixed and 𝑠 and 𝑏 1 may vary. By setting the right-hand sides of Eqns. ( 1)-( 3) to 0, we find that FP is determined by the only real root 𝑥 FP of the polynomial equation

𝑠 𝑘 𝑎 𝑥 3 -𝑘 (𝑠 + 1) 𝑥 2 + 𝑠 𝑎 1 𝑏 𝑥 + 𝑏 1 𝑏 = 0. (13)
The other coordinates of FP are given by

𝑦 FP = 𝑥 2 FP , and
𝑧 FP = 𝑠 𝑎 𝑥 3 FP -𝑥 2 FP (𝑠 + 1) 𝑏 = 𝑠 𝑎 1 𝑥 FP + 𝑏 1 𝑘 .
Note that for fixed values of 𝑎, 𝑏, 𝑎 1 , and 𝑘 the fixed point FP depends on 𝑠 and 𝑏 1 but not on 𝜖. However, its stability does depend on 𝜖 [START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF]. The Jacobian matrix of Eqns. ( 1)-( 3) is given by

J = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂𝑓 ∂𝑥 ∂𝑓 ∂𝑦 ∂𝑓 ∂𝑧 ∂𝑔 ∂𝑥 ∂𝑔 ∂𝑦 0 𝜖 ∂ℎ ∂𝑥 0 𝜖 ∂ℎ ∂𝑧 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
The stability of FP is determined by the eigenvalues of J evaluated at FP. Hence, the eigenvalues are the roots of the characteristic polynomial

P(𝜆) = ( 𝜆 -𝜖 ∂ℎ ∂𝑧 ) det ⎛ ⎜ ⎜ ⎝ 𝜆 - ∂𝑓 ∂𝑥 - ∂𝑓 ∂𝑦 - ∂𝑔 ∂𝑥 𝜆 - ∂𝑔 ∂𝑦 ⎞ ⎟ ⎟ ⎠ - 𝜖 ∂ℎ ∂𝑥 det ⎛ ⎜ ⎜ ⎝ - ∂𝑓 ∂𝑦 - ∂𝑓 ∂𝑧 𝜆 - ∂𝑔 ∂𝑦 0 ⎞ ⎟ ⎟ ⎠ ⇔ P(𝜆) = ( 𝜆 -𝜖 ∂ℎ ∂𝑧 ) det (𝜆 I -J f ) -𝜖 ∂ℎ ∂𝑥 ∂𝑓 ∂𝑧 ( 𝜆 - ∂𝑔 ∂𝑦 ) , (14)
where J f is the Jacobian matrix of the fast subsystem Eqns.

We are primarily interested in the stability of FP for 0 < 𝜖 ≪ 1 when FP is located on the unstable segments of the Z-shaped equilibrium curve of the fast subsystem, that is, on the saddle segment with 𝜆 1 < 0 < 𝜆 2 that lies between the two knees marked by saddle-node bifurcations (labeled SN 1 and SN 2 in Fig. 2) or on the segment of the upper branch with 𝜆 1 , 𝜆 2 > 0 that lies after the Hopf bifurcation (labeled HB 1 in Fig. 2). Note that the full system Eqns. ( 1)-( 3) is degenerate at 𝜖 = 0, because the Jacobian matrix J is then always singular. In the limit 𝜖 = 0 Eqn. ( 14) reduces to P(𝜆) = 𝜆 det (𝜆 I -J f ) , so that the eigenvalues of FP converge to the two eigenvalues 𝜆 1,2 of J f and the eigenvalue 𝜆 3 = 0 as 𝜖 → 0. At 𝜖 = 0, the characteristic polynomial P(𝜆) has two zero eigenvalues when FP lies on SN 1 (also 𝜆 1 = 0) or SN 2 (also 𝜆 2 = 0), and a pair of purely imaginary eigenvalues 𝜆 1 = λ2 along with 𝜆 3 = 0 when FP lies on HB 1 .

If FP lies on the lower branch, before SN 2 , then it is stable for 0 < 𝜖 ≪ 1, that is, 𝜆 3 < 0. As we decrease 𝑏 1 so that FP moves around the knee past SN 2 , the eigenvalues 𝜆 2 < 0 and 𝜆 3 < 0 become complex conjugate and move through the imaginary axis, after which they become real again and 𝜆 2 , 𝜆 3 > 0 (41). Hence, a Hopf bifurcation occurs that involves eigenvalues of both the fast (𝜆 2 ) and the slow (𝜆 3 ) equations. This is a singular Hopf bifurcation [START_REF] Braaksma | Singular Hopf bifurcation in systems with fast and slow variables[END_REF][START_REF] Guckenheimer | Singular Hopf bifurcation in systems with two slow variables[END_REF]. The characteristic polynomial Eqn. [START_REF] Chay | On the effect of the intracellular calcium-sensitive K + channel in the bursting pancreatic 𝛽-cell[END_REF] for FP at SN 2 , which corresponds to 𝑏 1 = 0 and FP = (0, 0, 0) in Fig. 2 for both 𝑠 = -1.61 and 𝑠 = -2.6, becomes P(𝜆) = (𝜆 + 1) (𝜆 2 + 0.2 𝜖 𝜆 + 𝜖 𝑎 1 𝑠), and the eigenvalues are 𝜆 1 = -1 and 𝜆 2,3 = -0.1 𝜖 ± 0.1 √ 𝜖 2 -100 𝜖 𝑎 1 𝑠. Therefore, provided 𝜖 > 0 is sufficiently small, 𝜆 2,3 are complex conjugate with negative real parts. Thus, the singular Hopf bifurcation does not happen exactly when FP lies on SN 2 , but when FP moves slightly past SN 2 onto the saddle-segment of the Z-shaped equilibrium curve. Indeed, the theory predicts that the singular Hopf bifurcation lies 𝑂(𝜖) close to the point where FP passes through the knee [START_REF] Guckenheimer | Singular Hopf bifurcation in systems with two slow variables[END_REF].

As we continue to decrease 𝑏 1 and trace FP along the saddle-segment of the Z-shaped equilibrium curve, the eigenvalues of FP are 𝜆 1 < 0 < 𝜆 2 and 𝜆 3 > 0 (Eqn. ( 15)). Note that as FP reaches the homoclinic bifurcation point HC the sum

A c c e p t e d m a n u s c r i p t

𝜆 1 + 𝜆 2 of the eigenvalues that converge to those of the fast subsystem for 𝜖 → 0 is negative for 𝑠 = -1.61 and positive for 𝑠 = -2.6. This is consistent with the type of the Hopf bifurcation in the fast subsystem, which is supercritical in the case of fold-homoclinic bursting (𝑠 = -1.61) and subcritical in the case of fold-subHopf (𝑠 = -2.6).

When FP moves around the right knee past SN 1 , one typically expects another singular Hopf bifurcation [START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model: A case study. A tutorial[END_REF]. However, we fixed 𝑎 = 0.5, 𝑏 = 1, 𝑎 1 = -0.1, and 𝑘 = 0.2 such that there is only one Hopf bifurcation on the upper branch of the Z-shaped equilibrium curve. This means that the equilibria of the fast subsystem are unstable on both sides of SN 1 . For 𝑏 1 such that FP is located 𝑂(𝜖) away from SN 1 on the middle branch of the Z-shaped curve, the eigenvalues of FP are 𝜆 1 < 0 < 𝜆 2 and 𝜆 3 > 0 (Eqn. ( 15)), while for FP located 𝑂(𝜖) away from SN 1 on the upper branch of the Z-shaped curve, its eigenvalues are 𝜆 1 , 𝜆 2 > 0 and 𝜆 3 < 0 (Eqn. ( 15)). Hence, in our case, it is not possible to get a singular Hopf bifurcation. Numerical calculations for 𝑠 = -2.6 and 𝜖 = 10 -6 seem to indicate that 𝜆 1 and 𝜆 3 pass through zero simultaneously, though generically, one would expect that they do so via two subsequent saddle-node bifurcations. However, none of the eigenvalues ever become zero (det J ∕ = 0) when 𝜖 > 0 and either 𝑠 = -2.6 or 𝑠 = -1.61; indeed, the degeneracy at SN 1 needs further investigation, but this is beyond the scope of this paper.

When 𝑏 1 is such that FP lies between HB 1 and SN 1 then FP is a saddle with two unstable eigenvalues (𝜆 3 < 0 for 0 < 𝜖 ≪ 1 Eqn. ( 15)). As 𝑏 1 decreases from a value with FP located close to SN 1 , the eigenvalues 𝜆 1 and 𝜆 2 that correspond to the fast subsystem coalesce on the real axis and become complex conjugate with positive real parts. This marks a transition from saddle-node to saddle-focus equilibrium. Finally, another Hopf bifurcation occurs as FP passes through HB 1 , but this is not a singular Hopf bifurcation because it involves only eigenvalues corresponding to the fast subsystem.

A.3 Computational method

The models simulations were done with the software package XPPAUT [START_REF] Ermentrout | Simulating, Analyzing and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students[END_REF]. The bifurcation analysis was performed with AUTO [START_REF] Doedel | AUTO97: continuation and bifurcation software for ordinary differential equations[END_REF].

  ). Published experimental data showing simultaneous recordings of voltage (𝑉 m ) and cytosolic calcium concen-2+ ] i ) in pancreatic 𝛽-cells (see Fig.2in (5)) and pituitary cells (see Fig.5in[START_REF] Van Goor | Dependence of pituitary hormone secretion on the pattern of spontaneous voltage-gated calcium influx. cell type-specific action potential secretion coupling[END_REF] and Fig.1in (47)) do, in fact, show that calcium levels in pancreatic 𝛽-cells oscillate between 0.15 𝜇M and 0.35 𝜇M, whereas in pituitary cells [Ca 2+ ] i can exceed 1 𝜇M during bursts.

  (1)-(2) for 𝜖 = 0 are plotted A c c e p t e d m a n u s c r i p t as well; for the logarithmic-scale pictures in Figs. 5(a) and 6(a) they are projected onto the (𝑧, 𝑥)-plane at the arbitrary values of 𝜖 = 10 -4 .

Figure 10

 10 point HB 2 in dependence on 𝑏 1 and 𝜖. The curves for both 𝑠 = -1.61 and 𝑠 = -2.6 are shown.

A c c e p t e d m a n u s c r i p t determined by the

  point where the 𝑧-nullcline intersects the Z-shaped equilibrium curve of the fast subsystem Eqns. (

A c c e p t e d m a n u s c r i p t

  (1)-(2) evaluated at FP. For fixed 𝑎 = 0.5, 𝑏 = 1, 𝑎 1 = -0.1, and 𝑘 = 0.2 the characteristic polynomial Eqn. (14) is given byP(𝜆) = (𝜆 + 0.2 𝜖) [ (𝜆 + 1) (𝜆 + 2 𝑠 𝑥 FP -1.5 𝑠 𝑥 2 FP ) + 2 𝑥 FP ]+𝜖 𝑎 1 𝑠 (𝜆 + 1).

Figure 1 .

 1 Figure 1. Bifurcation diagrams of the fast subsystem in the generic endocrine model (Appendix A.1) showing the bifurcations associated with square-wave (fold-homoclinic) bursting (panel a) and pseudo-plateau (fold-subHopf) bursting (panel b). (a) 𝑉 m L = -22.5 mV; (b) 𝑉 m L = -27.5 mV; HB -Hopf bifurcation; SN -saddle-node bifurcation; SNP -saddle-node of periodics; HC -homoclinic bifurcation point. Dashed lines denote instability. Sample bursting trajectories are superimposed on the bifurcation diagrams for each of the models.

Figure 2 .

 2 Figure 2. Bifurcation diagrams of the fast subsystem in the polynomial model Eqns. (1)-(2) showing the bifurcations associated with the transition from square-wave (fold-homoclinic) to pseudo-plateau (fold-subHopf) bursting as well as the 𝑧-nullclines (diagonal purple, blue and green lines) for three different values of the parameter 𝑏 1 and with (a) 𝑠 = -1.61, or (b) 𝑠 = -2.6; HB -Hopf bifurcation; SN -saddle-node bifurcation; HC -homoclinic bifurcation point. Dashed lines denote instability.

Figure 3 .

 3 Figure 3. (a) Three-dimensional view (𝜖, 𝑧, 𝑥) of the one-parameter bifurcation diagram with respect to 𝜖 of the full polynomial model Eqns. (1)-(3) in the case of square-wave bursting (𝑠 = -1.61, 𝑏 1 = -0.01); HB -Hopf bifurcation; SNP -saddle-node of periodics; FP -fixed point; HC -homoclinic bifurcation point. Dashed lines denote instability; (b) Sample bursting trajectories with increasing number of spikes, i.e., decreasing values of 𝜖 = 0.009; 0.005; 0.004; 0.0035; 0.0027; 0.0024; 0.002 are superimposed on the bifurcation diagram.

Figure 4 .

 4 Figure 4. (a) Bifurcation diagram of the full polynomial model Eqns. (1)-(3) with respect to 𝜖 in the case of pseudo-plateau bursting (𝑠 = -2.6, 𝑏 1 = -0.01); HB -Hopf bifurcation; SNP -saddle-node of periodics; FP -fixed point; PD -period-doubling bifurcation; HC -homoclinic bifurcation point. Dashed lines denote instability; (b) Sample bursting trajectories with increasing number of spikes, i.e., decreasing values of 𝜖 = 0.08; 0.06; 0.035; 0.023 are superimposed on the bifurcation diagram.

Figure 5 .A c c e p t e d m a n u s c r i p t 𝜖

 5t Figure 5. (a) Bifurcation of the full polynomial model Eqns. (1)-(3) in the case of square-wave bursting (𝑠 = -1.61, 𝑏 1 = -0.045); HB -Hopf bifurcation; PD -period-doubling bifurcation; FP -fixed point; HC -homoclinic bifurcation point. Dashed lines denote instability; (b) Sample bursting trajectories with increasing number of spikes, i.e., decreasing values of

Figure 6 .

 6 Figure 6. (a) Bifurcation diagram of the full polynomial model Eqns. (1)-(3) in the case of pseudo-plateau bursting (𝑠 = -2.6, 𝑏 1 = -0.21); HB -Hopf bifurcation; PD -period-doubling bifurcation; FP -fixed point; SNPsaddle-node of periodics; HC -homoclinic bifurcation point. Dashed lines denote instability; (b) Sample bursting trajectories with increasing number of spikes, i.e., decreasing values of 𝜖 = 1.105; 0.9; 0.6; 0.4; 0.25; 0.15 are superimposed on the bifurcation diagram.

Figure 7 .

 7 Figure 7. Simulations of the polynomial model showing the apparently chaotic spike-adding transition in the pseudo-plateau bursting regime in the case when (a) FP is well below the HC (𝑠 = -2.6, 𝑏 1 = -0.01); and (b) when FP is well above the HC (𝑠 = -2.6, 𝑏 1 = -0.21).

Figure 8 .

 8 Figure 8. Bifurcation diagrams of the full polynomial model Eqns. (1)-(3) in the cases of (a) square-wave bursting (𝑠 = -1.61, 𝑏 1 = -0.024), where FP lies just below the HC; and (b) pseudo-plateau bursting (𝑠 = -2.6, 𝑏 1 = -0.066), where FP lies slightly above the HC; HB -Hopf bifurcation; PD -period-doubling bifurcation; FP -fixed point; SNP -saddle-node of periodics; HC -homoclinic bifurcation point. Dashed lines denote instability.

Figure 9 .

 9 Figure 9. Simulations showing the behavior in the pseudo-plateau bursting regime for small 𝜖 in the cases (a) when FP is well below the HC (𝑠 = -2.6, 𝑏 1 = -0.01, 𝜖 = 0.0001); and (b) when FP is above the HC but below SN 1 (𝑠 = -2.6, 𝑏 1 = -0.12, 𝜖 = 0.001).

Figure 10 .

 10 Figure 10. Three-dimensional view (𝜖, 𝑧, -𝑏 1 ) of the two-parameter bifurcation diagram of the full polynomial model with respect to 𝜖 and 𝑏 1 showing the loci of HB 2 for 𝑠 = -1.61 (blue) and 𝑠 = -2.6 (red); HB -Hopf bifurcation; SN -saddle-node bifurcation. The vertical, dashed lines show the 𝑧-values of the hopf bifurcation HB 1 of the fast subsystem (𝜖 = 0).

Figure 2 :

 2 Figure 2:

Figure 3 :

 3 Figure 3:

Figure 4 :

 4 Figure 4:

Figure 5 :

 5 Figure 5:

Figure 6 :

 6 Figure 6:

Figure 8 :

 8 Figure 8:

2 HB 1 (Figure 10 :

 2110 Figure 9:

A c c e p t e d m a n u s c r i p t

  𝑠 𝑎 1 𝑥 + 𝑏 1 -𝑘 𝑧 is linear in 𝑥 and 𝑧 and represents Ca 2+ dynamics, with the term 𝑠 𝑎 1 𝑥 + 𝑏 1 replacing the source of calcium through voltage-gated calcium Ca 2+ channels and -𝑘 𝑧 standing for the decay term in Eqn.[START_REF] Belan | Calcium clamp in isolated neurones of the snail Helix pomatia[END_REF]. C4 The time-scale paramters 𝜙 and 𝜖 are such that 𝑥 and 𝑦 vary on a faster time scale than 𝑧. Although strictly speaking in the biophysical system Eqns. (4)-(6) there are three different intrinsic times scales for 𝑉 m , 𝑛 and 𝑐, it has usually been assumed that 𝑉 m and 𝑛 are fast variables compared to 𝑐. This is a reasonable assumption given that the time scale of change in Ca 2+ concentrations is several orders of magnitude smaller than 𝑉 m and 𝑛. Therefore, in the polynomial model we take 𝜙 = 1 and consider 𝜖 a small positive parameter.

	outward voltage sensitive K + currents; and -𝑏 𝑧 stands for the contribution
	of the outward calcium-sensitive potassium current.
	C2 The function 𝑔(𝑥, 𝑦) = 𝑥 2 -𝑦 is a quadratic function that gives a parabolic
	𝑦-nullcline and replaces the delayed rectifier activation kinetics (Eqn. (5) in
	Appendix A.1). It depends only on the membrane potential and is, therefore,
	decoupled from the third Eqn. (3).
	C3 The function ℎ(𝑥, 𝑧) = C5 The parameters 𝑎, 𝑏 ≥ 0 in the fast subsystem Eqns. (1)-(2) of the polynomial
	model are chosen, without loss of generality (W.L.O.G.), such that for a range
	of values of 𝑧 ≥ 0 there are three equilibrium points (𝑥 𝑒 𝑖 , 𝑦 𝑒 𝑖 ), 𝑖 = 1, 2, 3,
	given by the points of intersection of the 𝑥-and 𝑦-nullclines. Furthermore,
	we require that these equilibria are of the following type: (𝑥 𝑒 1 , 𝑦 𝑒 1 ) is a focus, (𝑥 𝑒 2 , 𝑦 𝑒 2 ) is a saddle, and (𝑥 𝑒 3 , 𝑦 𝑒 3 ) is a stable node. These conditions ensure
	that the fast subsystem Eqns. (1)-(2) of the polynomial model has a Z-shaped equilibrium curve defined by {𝑦 = 𝑥 2 and 𝑧 = ( 𝑠 𝑎 𝑥 3 -(𝑠 + 1) 𝑥 2 )
	is a cubic function that
	guarantees an N-shaped 𝑥-nullcline. Since 𝑥(𝑡) acts in place of the membrane
	potential (Eqn. (4) in Appendix A.1), the term -𝑠 (-𝑎 𝑥 3 + 𝑥 2 ) represents the
	contribution of the Ca 2+ inward current; -𝑦 represents the contribution of the

TABLE A

 A

	.1
	Parameter values of the Generic Pituitary Model
	𝑘 PMCA	20 s -1	𝑓 c	0.01
	𝑑 cell	10 𝜇m	𝑔 K(Ca)	0.2 nS
	𝑉 K	-65 mV	𝑔 Ca	0.81 nS
	𝑉 Ca	0 mV	𝑔 K	2.25 nS
	𝑉 m L	-27.5 mV	𝑉 n	0 mV
	𝑠 𝑚	12 mV	𝑠 𝑛	8 mV
	𝑘 s	1.25 𝜇M	𝜏 n	0.03 s -1
	A.2

Linear Stability Analysis of the Polynomial Model

  𝑥 FP , 𝑦 FP , 𝑧 FP ) that exists for all 𝜖 ≥ 0. The location of FP is

	According to Condition C6, the polynomial model Eqns. (1)-(3) has a unique equi-
	librium FP = (

Acknowledgments. AS was supported by the Intramural Research Program of NIDDK, National Institutes of Health, USA. HMO was supported by an Advanced Research Fellowship of the Engineering and Physical Sciences Research Council (EPSRC), UK. TR was supported by EPSRC grant EP/C544048/1. Both HMO and TR also received support from Cornell University via an IGERT grant of the National Science Foundation.

A Models and Methods

A.1 Generic Endocrine Model

The equations for the generic Hodgkin-Huxley-type model are:

where 𝐶 m is the membrane capacitance; 𝜏 n is the activation time constant for the delayed rectifier channel; 𝑛 ∞ is the steady state function for the activation variable 𝑛; 𝐶 m = 10 -5 × 𝐴 cell is the membrane capacitance; 𝑓 c is the fraction of free to total cytosolic Ca 2+ ; 𝛼 = 10 5 (2 × 9.65 × 𝐴 cell ) -1 is a factor that converts