The dynamics underlying pseudo-plateau bursting in a pituitary cell model

Pituitary cells of the anterior pituitary gland secrete hormones in response to patterns of electrical activity. Several types of pituitary cells produce short bursts of electrical activity which are more effective than single spikes in evoking hormone release. These bursts, called pseudo-plateau bursts, are unlike bursts studied mathematically in neurons (plateau bursting) and the standard fast-slow analysis used for plateau bursting is of limited use. Using an alternative fast-slow analysis, with one fast and two slow variables, we show that pseudo-plateau bursting is a canard-induced mixed mode oscillation. Using this technique, it is possible to determine the region of parameter space where bursting occurs as well as salient properties of the burst such as the number of spikes in the burst. The information gained from this one-fast/two-slow decomposition complements the information obtained from a two-fast/one-slow decomposition.


Introduction
Bursting is a common pattern of electrical activity in excitable cells such as neurons and many endocrine cells. Bursting oscillations are characterized by the alternation between periods of fast spiking (the active phase) and quiescent periods (the silent phase), and accompanied by slow variations in one or more slowly changing variables, such as the intracellular calcium concentration. Bursts are often more efficient than periodic spiking in evoking the release of neurotransmitter or hormone [1][2][3].
The endocrine cells of the anterior pituitary gland display bursting patterns with small spikes arising from a depolarized voltage [2][3][4][5]. Similar patterns have been observed in single pancreatic β-cells isolated from islets [6][7][8]. Figure 1(a) shows a representative example from a GH4 pituitary cell. Several mathematical models have been developed for this bursting type [5,[8][9][10]. Prior analysis showed that the dynamic mechanism for this type of bursting, called pseudo-plateau bursting, is significantly different from that of square-wave bursting (also called plateau bursting) which is common in neurons [11][12][13]. Yet this analysis did not determine the possible number of spikes that occur during the active phase of the burst. The goal of this paper is to understand the dynamics underlying pseudo-plateau bursting, with a focus on the origin of the spikes that occur during the active phase of the oscillation.
Minimal models for pseudo-plateau bursting can be written as 1V = f (V , n, c) (1.1) where V is the membrane potential, n is the fraction of activated delayed rectifier K + channels, and c is the cytosolic free Ca 2+ concentration. The velocity functions are nonlinear, and 1 and 2 are parameters that may be small. The variables V , n and c vary on different time scales (for details, see Section 2). By taking advantage of time-scale separation, the system can be divided into fast and slow subsystems. In the standard fast/slow analysis one considers 2 ≈ 0, so that V and n form the fast subsystem and c represents the slow subsystem. One then studies the dynamics of the fast subsystem with the slow variable treated as a slowly varying parameter [12,[15][16][17][18]. This approach has been very successful for understanding plateau bursting, such as occurs in pancreatic islets [19], pre-Bötzinger neurons of the brain stem [20], trigeminal motoneurons [21] or neonatal CA3 hippocampal principal neurons [14], Fig. 1(b). It has also been useful in understanding aspects of pseudo-plateau bursting such as resetting properties [11], how fast subsystem manifolds affect burst termination [17], and how parameter changes convert the system from plateau to pseudo-plateau bursting [12].
An alternate approach, which we use here, is to consider 1 ≈ 0, so that V is the sole fast variable and n and c form the slow subsystem. With this approach, we  [14].
show that the active phase of spiking arises naturally through a canard mechanism, due to the existence of a folded node singularity [22][23][24][25]. Also, the transition from continuous spiking to bursting is easily explained, as is the change in the number of spikes per burst with variation of conductance parameters. Thus, the one-fast/twoslow variable analysis provides information that is not available from the standard two-fast/one-slow variable analysis in the case of pseudo-plateau bursting.

The mathematical model
We use a model of the pituitary lactotroph, which produces pseudo-plateau bursting over a range of parameter values [10]. To achieve a minimal form, we use the model without A-type K + current (I A ). It includes three variables: V (membrane potential), n (fraction of activated delayed rectifier K + channels), and c (cytosolic free Ca 2+ concentration). The equations are: where I Ca is an inward Ca 2+ current, I K is an outward delayed rectifying K + current, I K(Ca) is a small-conductance Ca 2+ -activated K + current, and I BK is a fastactivating large-conductance BK-type K + current. The currents in the equations above are: (2.7) The steady state activation functions are given by: Default parameter values are given in Table 1. The variables V , n and c vary on different time scales. The time constant of V is given by . During a bursting oscillation, the minimum of g T otal is 0.483 pS and the maximum is 3 pS. Hence, min g T otal , or 1.7 ms ≤ τ V ≤ 10.4 ms, for C m = 5 pF, a typical capacitance value for lactotrophs. The time constant for n is τ n = 43 ms. For the variable c, the time constant is 1 f c k c = 1 (0.01)(0.16) ms = 625 ms. Thus, n and c change more slowly than V . This time scale separation between V and (c, n) can be accentuated when C m is made smaller than the default 5 pF, i.e., when C m → 0, τ V gets smaller and V varies much faster. Thus, we can view the capacitance C m as a representative of the dimensionless singular perturbation parameter 1 in this model (Eq. 1.1).
All numerical simulations and bifurcation diagrams (both one-and two-parameter) are constructed using the XPPAUT software package [26], using the Runge-Kutta integration method, and computer codes can be downloaded from the following website: http://www.math.fsu.edu/~bertram/software/pituitary. The surface in Fig. 9 was constructed using the AUTO software package [27]. All graphics were produced with the software package MATLAB.

The reduced system
We consider the full system (Eqs. (2.1)-(2.3)) as having one fast variable V and two slower variables n and c. The time-scale separation can be accentuated by decreasing the singular perturbation parameter C m . This facilitates analysis of the system dynamics [28]. In the limit C m → 0, the trajectories of the system lie on a 2-D surface called the critical manifold. If we define the right hand side of Eq. (2.1) by then the critical manifold is the surface S satisfying 2) The equation f (V , c, n) = 0 can be solved in explicit form for n as The critical manifold (3.3) is a folded surface ( Fig. 2) that consists of three sheets separated by two fold curves (L − and L + ). The lower and upper sheets are attracting ( ∂f ∂V < 0) and the middle sheet is repelling ( ∂f ∂V > 0). The lower (L − ) and upper (L + ) fold curves are given by This yields two constant V values and two equations for n in the form of n = n(c). Thus, the fold curves (L ± ) are (V ± , c, n ± (c)) where V − and V + are constant V values. The curve L + is projected vertically (along the fast variable V ) onto the lower sheet to obtain the projection curve P(L + ), and similarly for the (L − ) projection onto the upper sheet. Figure 2 shows the critical manifold, the fold curves and the projections of the fold curves. The reduced flow (when C m → 0) is described by (3.3), the differential equation for c (Eq. (2.3)), and a differential equation for V which can be obtained by differentiating f (V , c, n) = 0 with respect to time. That is, where n satisfies Eq. (3.3), andṅ,ċ satisfy Eqs. (2.2), (2.3). The two differential equations for the reduced system are thus Since ∂f ∂V = 0 on L ± , the reduced system is singular along the fold curves. The system can be desingularized by rescaling time with τ = −( ∂f ∂V ) −1 t. The desingularized system is then we have the desingularized system (3.12) The desingularized system describes the flow on the critical manifold. Because of the time rescaling, the flow on the middle sheet, where ∂f ∂V > 0, must be reversed to obtain the equivalent reduced flow.

Folded singularities and canards
Equilibria of the desingularized system are classified as ordinary singularities and folded singularities. An ordinary singularity is an equilibrium of Eqs. (2.1)-(2.3) and satisfies A folded singularity lies on a fold curve (L + or L − ), and satisfies: A folded singularity is classified as a folded node if the eigenvalues are real and have the same sign, a folded saddle if the eigenvalues are real and have opposite signs, or a folded focus if the eigenvalues are complex [22,23,25,29]. For parameter values used in Fig. 2, the system has a folded node (with negative eigenvalues) on L + (FN, blue point, in Fig. 2), and a folded focus on L − (not shown).
There are an infinite number of singular trajectories on the top sheet that pass through the folded node (FN). These are called singular canards [22]. The singular canard that enters the FN in the direction of the strong eigenvector is called the strong canard (SC, green curve, in Fig. 2). This curve and the fold curve L + delimit the singular funnel that consists of all initial conditions whose trajectories for the reduced system pass through the folded node. The singular funnel and key curves are projected onto the (c, V )-plane in Fig. 3. The different panels are obtained with different values of the parameter g K .

Singular periodic orbits, relaxation oscillations, and mixed mode oscillations
A singular periodic orbit (Fig. 2, black curve with arrows) can be constructed by solving the desingularized system for the flow on the top and bottom sheets of the critical manifold, and then projecting the trajectory from one sheet to the other along fast fibers when the trajectory reaches a fold curve. The singular periodic orbit is the closed curve constructed in this way. This process was discussed in detail in [22,28,30]. Briefly, the trajectory moves along the bottom sheet until L − is reached. At this point the reduced flow is singular ( ∂f ∂V = 0). The quasi-steady state assumption f (V , c, n) = 0 is no longer valid and there is a rapid motion away from the fold curve L − . This rapid motion is seen as vertical movement to the top sheet (the dynamics are governed by the layer problem, see [22,28]). The trajectory moves to a point on P(L − ) and from there is once again governed by the desingularized equations, moving along the top sheet until L + is reached. The fast vertical downward motion along fast fibers returns the trajectory to a point on P(L + ) on the bottom sheet, completing the cycle.
When the singular periodic orbit reaches L − it jumps up to a point on P(L − ). If this point on P(L − ) is in the singular funnel, then the orbit will move through the FN. Otherwise it will not. Let δ denote the distance measured along P(L − ) from the phase point on P(L − ) of the singular periodic orbit to the strong canard (SC in Fig. 3). When the phase point is on the strong canard, δ = 0. Let δ > 0 when the phase point is in the singular funnel and δ < 0 when the phase point is outside the singular funnel. Singular canards are produced when δ > 0.
In Fig. 3(a) the singular periodic orbit jumps to a point on P(L − ) outside of the singular funnel (δ < 0), so it does not enter the FN. This orbit is a relaxation oscillation [31]. In Fig. 3(b) δ > 0, so the orbit is a singular canard. Away from the singular limit, this singular canard perturbs to an actual canard that is characterized by small oscillations about L + [22]. The combination of these small oscillations with the large oscillations that occur due to jumps between upper and lower sheets yields mixed mode oscillations [24,32]. The small oscillations have zero amplitude in the singular case, which grows as √ C m for C m sufficiently small [23]. A discriminating condition between relaxation and mixed mode oscillations is δ = 0, where the singular periodic orbit jumps to P(L − ) on the SC curve.
When C m > 0 the full system (Eqs. (2.1)-(2.3)) produces spiking for δ < 0 and mixed mode oscillations for δ > 0. Figure 4 shows these two different cases for g BK = 0.4 nS. For g K = 5.1 nS (δ < 0 in Fig. 3(a)), the nearly-singular periodic orbit produced when C m = 0.001 pF ( Fig. 4(a)) perturbs to continuous spiking when C m = 10 pF (Fig. 4(e)). When g K = 4 nS the singular periodic orbit enters the singular funnel (Fig. 3(b)), so when C m is increased the singular orbit transforms to mixed mode oscillations. For C m = 0.5 pF mixed mode oscillations with small spikes are produced ( Fig. 4(d)). As C m is increased to 10 pF, mixed mode oscillations with larger spikes are produced. This is the genesis of pseudo-plateau bursting ( Fig. 4(f)).

Analysis of the desingularized system and folded nodes
We next discuss the singularities of the desingularized system for a range of g K and g BK values (Fig. 5). The system (with g BK = 0.4 nS) has a single-branched Vnullcline (green curve) that satisfies F (V , c, n) = 0 and a three-branched c-nullcline (orange curves) L − , L + and CN1. The curves L − , L + satisfy ∂f ∂V = 0, and are the same as the fold curves in Fig. 3. The curve CN1 satisfies αI Ca + k c = 0. There are folded singularities that are located at intersections of the V -nullcline with L − or L + , and ordinary singularities that are located at intersections with CN1. For fixed g BK , changing g K affects the position of the V -nullcline but not the c-nullcline.
For values g K < 0.5131 nS, there is a stable node on CN1 (A 1 ), which would be on the top sheet of the critical manifold. There are also two folded saddles on L + (B 1 . Filled circles represent stable singularities and unfilled circles represent unstable singularities. Red circles (filled or unfilled) are ordinary singularities. Filled and unfilled circles in blue are folded nodes and folded saddles, respectively. Filled circles in cyan are folded foci. The points TR1 and TR2 are transcritical bifurcations (type II folded saddle-node bifurcations) and SN1 and SN2 are standard saddle-node bifurcations (type I folded saddle-node bifurcations). and C 1 ) and two folded foci on L − (D 1 and E 1 ). When g K is increased to 0.5131 nS the stable node A 1 moves down and to the left and the folded saddle B 1 moves to the left. These two equilibria coalesce at a transcritical bifurcation (TR1). This transcritical bifurcation corresponds to a bifurcation of folded singularities called a type II folded saddle-node [22,30,33]. Following this bifurcation, the folded singularity is a folded node. For g K = 4 nS, the equilibria on L + are the folded node (B 3 ) and the folded saddle (C 3 ). The equilibrium on CN1 (A 3 ) is now a saddle point. There is no qualitative change of equilibria on L − . When g K is increased to 7.588 nS the equilibria B 3 and C 3 coalesce at a saddlenode bifurcation point (SN1). This is a standard saddle-node bifurcation of folded singularities and is called a type I folded saddle-node [22,30,33]. As g K is increased to 43.1 nS, the folded focus D 5 moves to the left and changes to a folded node at D 6 . The saddle points on CN1 move downward and to the left as g K is increased. For g K = 129.2 nS, the saddle point A 6 coalesces with the fold node D 6 at a second transcritical bifurcation (TR2); again a type II folded saddle-node. Beyond this, the ordinary singularity (A 8 , A 9 ) is stable and the folded singularity becomes a folded saddle. Moreover the folded focus E 6 has become a folded node (E 7 ). As g K is in-creased further to 137.2 nS, there is a second type I saddle-node bifurcation (SN2) at which the folded node and the folded saddle coalesce and disappear. For the values g K > 137.2 nS, the only equilibrium is on CN1 and is an ordinary stable node (A 9 ). This is on the bottom sheet of the critical manifold.
Varying g BK slightly affects the V -nullcline and strongly affects the c-nullcline in the (c, V )-phase plane. Increasing g BK moves the fold curves together, eventually taking the fold out of the critical manifold. Figure 6 shows qualitative changes in the equilibria when g BK is varied, with g K = 7.588 nS. When g BK = 0.2 nS there is a saddle point on CN1 (A) and two folded foci (D and E) on L − (Fig. 6(a)). When g BK is increased to 0.4 nS, the curve L + moves down and a type I folded saddle-node bifurcation occurs (SN1 in Fig. 6(b)). When g BK is increased further, the saddle-node splits into a folded node (B) and a folded saddle (C) on L + , as shown for g BK = 1 nS in Fig. 6(c).
The folded node (B) and the saddle point (A) coalesce at a transcritical bifurcation (type II folded saddle-node) when g BK = 3.96 nS (TR1 in Fig. 6(d)). Beyond this, the ordinary singularity (A) is a stable node that lies on the top sheet of the critical manifold. When g BK = 20 nS the folded singularities are either saddles or foci, Fig. 6(e). For g BK ≈ 32.12 nS the two folded foci on L − change to folded nodes. Finally, when g BK is increased to 32.1224 nS, the fold curves L + and L − merge. As a result, the folded saddles coalesce with the folded nodes at type I folded saddlenode bifurcations (SN3 and SN4 in Fig. 6(f)). Beyond this, there is only a stable node (A in Fig. 6(g)). The disappearance of the L + and L − curves correspond to the disappearance of the fold in the critical manifold.
The two-parameter bifurcation diagram in Fig. 7 summarizes the variations of the bifurcations in Fig. 5 and Fig. 6 over a range of g K and g BK values. The curves TR1 and TR2 correspond to the transcritical bifurcations (type II folded saddle-node bifurcations), and SN1-SN4 correspond to the saddle-node bifurcations (type I folded saddle-node bifurcations). At g BK = 32.1224 nS the L + and L − lines coalesce into a single line. This contains the SN3 and SN4 bifurcations, up until SN3 and SN4 coalesce at a codimension-2 bifurcation (for g K = 83.7122 nS). For large g K , the L + /L − line contains no folded singularities (dashed line).
For g K and g BK values in regions A, D and E there is only a stable node and the full system is in a depolarized (A) or hyperpolarized (D or E) steady state. In the left portion of region C there is a folded focus which becomes a folded node in the right portion of C. This family of folded singularities is on L − . In region D there is a folded node on L − for negative values of c. Region B consists of the folded nodes on L + , and it is the key region for the existence of mixed mode oscillations, since δ > 0 for much of this region (shown below).

Twisted slow manifolds and secondary canards
The folded nodes discussed above are important since they yield small oscillations (for C m > 0) in all trajectories entering the singular funnel. In this section we explain the genesis of those oscillations (for more details, see [22,23,28,32]).  Fig. 5.  Fig. 7 Two-parameter bifurcation structure for the desingularized system. The curves TR1 and TR2 represent the transcritical bifurcations (type II folded saddle-nodes). The curves SN1-SN4 represent saddle-node bifurcations (type I folded saddle-nodes). The horizontal line is where the fold curves L + and L − coalesce. A codimension-2 bifurcation occurs at the intersection of the SN curves.
Folded nodes or saddles are characterized by the ratio of their eigenvalues. Let λ 1 and λ 2 be the eigenvalues of the folded singularity on the fold curve L + such that In region A of Fig. 7, which consists of folded saddles, μ < 0. On the TR1 curve μ = 0 since λ 1 = 0. Folded nodes occur in region B, so μ > 0. For C m > 0, but small, a trajectory approaching a folded node will oscillate, due to twists in the attracting and repelling sheets of the slow manifold. The maximum number of oscillations is given by [23,32] S max = μ + 1 2μ , (5.2) which is the greatest integer less than or equal to μ+1 2μ . At a point in region B and close to the TR1 curve in Fig. 7, μ > 0 but small. Hence, S max is large. Similarly on SN1 μ = 0, so in region B and close to the SN1 curve μ > 0 and small, so S max is large. Between these curves μ increases and S max declines. This is shown in Fig. 8 for the case g BK = 0.4 pS. The small value of μ over the full range of g K values in (Fig. 8a) suggests the system is close to a folded saddle-node bifurcation, either type I (SN1) or type II (TR1).
The attracting sheets of the critical manifold (S a ) and the repelling middle sheet (S r ) come together at the fold curves L + and L − . For C m > 0, Fenichel theory [34] tells us that the critical manifold perturbs smoothly to invariant attracting (S a,C m ) and repelling (S r,C m ) manifolds away from L + and L − . However, the critical manifold is non-hyperbolic on L + and L − , and perturbs to twisted sheets near these curves to preserve uniqueness of solutions [23,35]. Figure 9 shows how the top attracting S + a,C m (blue) and middle repelling S r,C m (red) sheets of the slow manifold intersect and twist. The numerical method used to compute the slow manifolds was developed by Desroches et al. [36,37].
The primary weak canard corresponds to the weak eigendirection of the folded node. It is at the intersection of the invariant manifolds S + a,C m and S r,C m and serves  , blue surface) and middle repelling (S r,C m , red surface) sheets of the slow manifold are twisted around the blue dashed curve, which is the axis of rotation. The primary strong canard (SC, green) moves from the attracting to the repelling sheet without any rotations. The secondary canards ξ 1 (gray curve, one rotation), ξ 2 (purple curve, two rotations) and ξ 3 (gold curve, three rotations) flow from the attracting to repelling sheet with different numbers of rotations. A portion of the pseudo-plateau burst trajectory (PPB, black curve) is superimposed and has two small oscillations. The full system has unstable equilibrium (cyan, filled curcle). as their axis of rotation. All other canards twist about the primary weak canard; they follow S + a,C m as it twists and then follow S r,C m for a distance as it twists. The primary strong canard, which corresponds to the strong eigendirection of the folded node, moves along S + a,C m to S r,C m without any rotation (SC, green curve in Fig. 9). Other, secondary, canards rotate a number of times, depending on how close they are to the primary strong canard. A secondary canard that makes k small rotations in the vicinity of the folded node is called the k th secondary canard. Figure 9 shows the first (ξ 1 , gray), second (ξ 2 , purple) and third (ξ 3 , olive) secondary canards that make one, two and three rotations, respectively. For C m > 0, but small, there are S max − 1 secondary canards which divide the funnel region between the primary canards into S max subsectors [24]. The first subsector is bounded by the strong canard SC and the first secondary canard ξ 1 and trajectories entering here have one rotation. The second subsector is bounded by ξ 1 and ξ 2 and trajectories entering here have two rotations. The last subsector is bounded by the last secondary canard and the primary weak canard. The maximal rotation number is achieved in the last subsector; trajectories entering here have S max rotations [23,28,32]. Figure 9 also shows a portion of the pseudo-plateau burst trajectory (PPB, black curve) for C m = 2 pF. It enters the funnel region in the rotational subsector bounded by ξ 1 and ξ 2 , and hence, makes two full rotations and then leaves the repelling sheet as it moves towards the lower attracting manifold S − a,C m . These rotations are the small oscillations or "spikes" during the burst active phase. Figure 10(a) shows burst trajectories for three values of C m projected onto the (c, V )-plane. Also shown are L + , L − , the singular strong canard SC and the folded node of the desingularized system. Finally, the line along the weak eigendirection of the folded node is included (WED, pink curve). With C m = 0.001 pF the system is nearly singular and the "bursting" trajectory enters and leaves the folded node along the WED. The small oscillations near the folded node are too small to see. The region near the folded node is magnified in Fig. 10(b). With C m = 0.1 pF the burst trajectory again passes through the folded node along the WED, but now the small oscillations are visible in Fig. 10(b). The small oscillations of this burst trajectory first decrease and then increase in amplitude. This is often seen in mixed mode oscillations that are associated with a folded node singularity, in contrast to those associated with a singular Hopf bifurcation, where the amplitude of successive small oscillations increases [38]. Finally, with C m = 2 pF (the value used in Fig. 9) the small oscillations are prominent even in the larger vertical scale used in Fig. 10(a).

The boundaries of mixed mode oscillations
For a periodic mixed mode oscillation (i.e., pseudo-plateau bursting) solution to exist, there must be a folded node singularity and the periodic orbit must enter the funnel. In this section we construct curves in the two-parameter g K -g BK plane that form boundaries for the existence of mixed mode oscillations.
From Fig. 7 we know that folded node singularities only occur in regions B and C (and in region D for negative values of the Ca 2+ concentration). Those in region C occur on L − and the periodic orbit never enters the corresponding singular funnel. We therefore focus on region B. This region is highlighted in Fig. 11(a). Above the TR1 curve the system has a depolarized stable steady state. Below the SN1 curve the system spikes continuously. Between these curves a folded node singularity exists, and the requirement for periodic mixed mode oscillations is that δ > 0. That is, the singular orbit must enter the singular funnel. Thus, the final curve delimiting the MMOs region is δ = 0 (the set of g K and g BK values at which the singular periodic orbit intersects both the strong canard and the curve P(L − )), shown in green in Fig. 11. For parameter values between the δ = 0 and TR1 curves periodic mixed mode oscillations, i.e., pseudo-plateau bursting, exist and are stable. This critical region is magnified in Fig. 11(b). Figure 12 shows how the burst duration and the number of spikes in a burst vary over a range of g K and g BK values for C m = 5 pF. A similar map of parameter space was used previously in the analysis of a parabolic burster [39]. Two-parameter Fig. 12 The active phase duration and the number of spikes per burst of the full system for C m = 5 pF. The system displays steady state, spiking or bursting solutions. Steady state and spiking solutions are represented by black dots and small black circles, respectively. The bursting region is bounded by the supercritical Hopf bifurcation (HB, black) and the right branch of the period doubling (PD, green) curves. The bursting patterns in this region are represented by colored circles. The size of a circle represents the active phase duration, with larger circles corresponding to longer active phase durations. The color of a circle represents the number of spikes per burst. Cyan circles correspond to smaller number of spikes per burst (minimum of two spikes) and the largest dark red circle corresponds to the largest number of spike per burst (36 spikes in a burst). bifurcation curves of the full system (Eqs. (2.1)-(2.3)) are also shown. These include a curve of supercritical Hopf bifurcations (HB, black) and a curve of period doublings (PD, green). To the left of the HB curve the system is at a steady state (black dots), and to the right of and above the PD curve the system produces continues spiking (small black circles). For the values of g K and g BK inside the PD curve the system produces pseudo-plateau bursting oscillations (MMOs), represented by colored circles.
In the bursting region the active phase duration and the number of spike per burst vary with respect to the values of g K and g BK . The size of each circle represents the active phase duration, and the color of the circle (from cyan to dark red) represents the number of spikes in a burst. A burst with larger number of spikes has longer active phase duration, and in an actual cell this determines the amount of Ca 2+ influx and hormone released. The bursts that have the shorter active phase duration and the smaller number of spikes occur near the right branch of the PD curve. These bursts are represented by smaller cyan circles in Fig. 12. For example, when g BK = 1 nS and g K = 6 nS the system produces bursting oscillations with three spikes per burst (as in Fig. 4(f)). When one moves away from the right to the left branch of the PD curve by increasing g BK or decreasing g K the burst duration becomes longer and the number of spikes in a burst becomes larger. The longest active phase duration is about 8.4 sec and the largest number of spikes per burst is about 36, represented by the largest dark red circle. These values will change when C m is changed.
The region between the HB and the left branch of the PD curves is bistable between bursting and continuous spiking. Orange circles with small black circles at the centers represent bistable solutions that are simulated by varying the initial condi-tions. This shows that the borders of the bursting region are the HB and the right branch of the PD curves. The dark blue circles represent bursting oscillations without small oscillations since the amplitudes of the spikes are almost zero, i.e., the small oscillations are too small to see.
The results that are shown in Fig. 12 are very consistent with the analysis of the mixed mode oscillations in Fig. 11. The HB and TR1 curves overlap, demonstrating that for small C m the HB of the full system corresponds to a type II saddle-node bifurcation of the desingularized system. Also, the HB curve and the left branch of the PD curve are almost indistinguishable for small C m . For these C m values (C m < 0.001 pF), the right branch of the PD curve converges to the δ = 0 curve of the desingularized system. Hence, the left and right borders of the MMOs in the singular limit C m → 0 pF correspond to the left and right borders of the bursting region of the full system for C m > 0, with the exception that the bursting region is smaller for larger values of C m . Also, the bistable region between the PD and HB curves only exists as the left PD moves away from the HB, which occurs as C m is increased.
In Fig. 11 the MMOs region delimited by the TR1 and δ = 0 curves can be divided into subregions that have different numbers of small oscillations. For parameter values in the subregion near the curve δ = 0 the periodic orbit enters the funnel region near the strong (primary) canard. This subregion corresponds to the first subsector of the funnel region, and for C m > 0 only one small oscillation occurs in a burst. This corresponds to the jump from the lower attracting sheet to the upper attracting sheet and is not due to the folded node. When one moves leftward by decreasing g K , δ increases and the periodic orbit enters the funnel region through other subsectors. As a result, the number of small oscillations in a burst increases. When one moves to the subregion near or on the TR1 curve by decreasing g K further, the periodic orbit enters the funnel region through the last subsector. The number of small oscillations is closer to S max , the maximum number of spikes in a burst as determined by the eigenvalues of the folded node. Moreover, increasing g BK has the same effect as decreasing g K . These trends in the number of small oscillations obtained from an analysis of the desingularized system [28] are expressed far from the singular limit as shown in Fig. 12 where C m = 5 pF. Here the longest bursts occur near the HB curves, as predicted.

A comparison with a two-fast/one-slow variable analysis
Using a one-fast/two-slow variable analysis we have shown the genesis of the spikes in a burst and how the number of spikes in a burst varies in the g K -g BK parameter space. The regions for steady states, pseudo-plateau bursting (mixed mode oscillations) and spiking are clearly identified in this parameter space (Fig. 11). This has been done by investigating the qualitative changes of the desingularized system when parameters g K (Fig. 5) and g BK (Fig. 6) are varied, which are summarized in Fig. 7.
Here we investigate whether this information can be obtained from a standard twofast/one-slow variable analysis. Figure 13(a) shows a bifurcation diagram of the V -n fast subsystem with c treated as a parameter (referred to as a "z-curve"). The subsystem is bistable over a large range of c values, with stable depolarized and hyperpolarized steady states, separated by saddle points. The c-nullcline is superimposed, now Two-fast/one-slow analysis for g BK = 0.4 nS, C m = 10 pF and different values of g K . The black "z-curve" is the curve of equilibria of the V -n fast subsystem. This has stable (solid) and unstable (dashed) branches. (a) g K = 0.1 nS, the full system with stable equilibrium (A 1 ) is in a depolarized steady state. (b) g K = 4 nS, the fast subsystem has an unstable limit cycle that emerges from the subcritical Hopf bifurcation (subHB). Pseudo-plateau bursting (PPB, black trajectory) is produced. The equilibrium of the full system (A 3 ) is unstable. (c) g K = 5.1 nS, the full system produces periodic spiking that appears unrelated to the fast subsystem bifurcation structure. The equilibrium point of the full system (A) is unstable. In all panels the system is bistable over a range of c-values. The points A 1 and A 3 are the same equilibrium points as A 1 and A 3 in Fig. 5, respectively. thinking of c as a slowly-changing variable rather than as a parameter. This is the standard approach used in a two-fast/one-slow variable analysis. In all three panels of Fig. 13 parameters are set at g BK = 0.4 nS, C m = 10 pF, and g K is varied.
In Fig. 13(a), with g K = 0.1 nS, there is an intersection of the c-nullcline on the upper stable branch at location A 1 . This is a stable equilibrium of the full 3dimensional system, and corresponds to A 1 in the analysis shown in Fig. 5. Thus, both types of analysis indicate that the system will come to rest at a depolarized steady state when g K = 0.1 nS.
When g K is increased there is a subcritical Hopf bifurcation on the upper branch with emergent unstable periodic solutions of the fast-subsystem. This is shown in Fig. 13(b) for the case g K = 4 nS. Pseudo-plateau bursting occurs for this and nearby values of g K . The full system unstable equilibrium (A 3 ) corresponds to A 3 in Fig. 5.
The superimposed burst trajectory in Fig. 13(b) only weakly follows the fastsubsystem bifurcation diagram. Most notably, there are no stable periodic solutions of the fast subsystem, only bistability between two steady states. Also, the trajectory never follows the lower branch of stationary solutions and greatly overshoots the lower knee.
The subcritical Hopf bifurcation migrates leftward when g K is increased to 5.1 nS. The unstable branch of periodics goes through a saddle-node bifurcation, yielding a branch of stable periodic solutions of the fast subsystem (Fig. 13(c)). There is bistability between upper and lower branches of the z-curve which is typically a necessary condition for bursting with this type of analysis. However, bursting is not produced for this value of g K . Instead, the system spikes continuously.
This example illustrates that features well described by the one-fast/two-slow variable analysis are not at all well described by a standard two-fast/one-slow variable analysis. Most notably, the transition from bursting to spiking is well characterized in the one-fast/two-slow variable analysis as the point at which δ = 0. Note that this is not a bifurcation point of the desingularized system, but reflects the jump point from the lower sheet of the slow manifold to the upper sheet. In contrast, the bursting to spiking transition is not predicted from the two-fast/one-slow analysis, and indeed the periodic spiking trajectory of the full system occurs over a range of the fast-subsystem bifurcation diagram that contains only stable equilibria. The one-fast/two-slow approximation is good even at higher values of C m , for example, when C m = 5 pF (Fig. 12). Similar remarks apply for smaller values of C m , where the one-fast/twoslow approximation becomes more accurate while the two-fast/one-slow approximation does not. The two-fast/one-slow approximation becomes more accurate when c is much slower than both V and n, but in this case only a stable steady solution or a relaxation oscillation is produced.

Discussion
The canard mechanism has been used to understand mixed mode oscillations in several neuronal models [30,37,[40][41][42][43][44]. In these examples, the small oscillations correspond to subthreshold oscillations that occur between the electrical impulses. We have previously analyzed pseudo-plateau bursting in a pituitary lactotroph model using canard theory [28]. However, the model used was a simplification in which the cytosolic free Ca 2+ concentration was treated as a fixed parameter and the second slow variable (in addition to the variable n used here) was an inactivation variable for an A-type K + current. In the current paper, we again focused on pseudo-plateau bursting in a pituitary lactotroph model, but now with emphasis on a BK-type K + current. In this analysis, we have examined the effects of changing the parameters C m , g K and g BK . The parameter g BK is important for producing bursting oscillations in actual pituitary cells in which bursting is converted to spiking when BK-type K + channels are blocked [45].
Here, using C m to control the separation in time scales, we identified two slow variables (n, c) and one fast variable (V ). Using the one-fast/two-slow variable analysis we showed that pseudo-plateau bursting is a canard-induced mixed mode oscillation. There are two main requirements for the existence of these bursting oscillations [22][23][24]32]. One is that the desingularized system must have a folded node singularity, i.e., the eigenvalue ratio (μ) has to be positive. The second requirement is that the singular periodic orbit should enter the singular funnel and pass through the folded node, i.e., δ should be positive. In short, canard-induced mixed mode oscillations exist if both μ and δ are positive.
Using this technique we can understand several features of the burst and several trends that occur as parameters are varied. When both μ and δ are positive, small oscillations are produced during the active phase of a burst and their amplitude is proportional to √ C m for C m sufficiently small [23]. We obtained the bursting borders in the (g K , g BK )-plane (Figs. 11 and 12), and predicted how the active phase duration and the number of spikes per burst vary with changes in parameters.
The singular perturbation analysis performed here is technically more effective and informative in the singular limit (i.e., for sufficiently small values of C m ) [22,23]. However, it provides useful information even far from this limit, as we showed in Figs. 11 and 12. Eventually, as the singular parameter (C m ) is increased sufficiently, new dynamics will be introduced, and the insights from the singular analysis are no longer valid.
The one-fast/two-slow decomposition used here contrasts with the two-fast/oneslow variable analysis used previously for pseudo-plateau bursting [10][11][12][13]. Our analysis explains the origin of the small-amplitude spikes that occur during the active phase of pseudo-plateau bursting, the transition between spiking and bursting, and information about how the number of spikes per burst varies with parameters. While the two-fast/one-slow variable analysis provides little information on these things, it does provide valuable information about how one can make a transition between plateau and pseudo-plateau bursting as one or more parameters are changed [12]. It also provides information about complex phase resetting properties [11] and the termination of spikes in a burst [17]. Both fast/slow decompositions are approximations, however, to a system that evolves on three time scales. Some studies [13,17,18] focus on the dynamics of the full system, and illustrate the complexity of the seemingly simple set of equations. The advantage of obtaining useful information of the full system by a two-fast/one-slow or one-fast/two-slow decomposition points to the fact that system (2.1)-(2.3) actually evolves on three time scales: V fast, n intermediate and c slow. This can also be seen by the magnitude of μ which is bounded from above by μ max ≈ 0.07 ( Fig. 8(a)). Hence, we are close to folded saddle-node regimes (type I and type II) [33,38] and a more detailed bifurcation analysis may explain the relation between the two-fast/one-slow and one-fast/two-slow splitting. This is left for future work.