Geometric slow–fast analysis of a hybrid pituitary cell model with stochastic ion channel dynamics

To obtain explicit understanding of the behavior of dynamical systems, geometrical methods and slow–fast analysis have proved to be highly useful. Such methods are standard for smooth dynamical systems and increasingly used for continuous, non-smooth dynamical systems. However, they are much less used for random dynamical systems, in particular for hybrid models with discrete, random dynamics. Here we pro-pose a geometrical method that works directly with the hybrid system. We illustrate our approach through an application to a hybrid pituitary cell model in which the stochastic dynamics of very few active large-conductance potassium (BK) channels is coupled to a deterministic model of the other ion channels and calcium dynamics. To employ our geometric approach, we exploit the slow–fast structure of the model. The random fast subsystem is analyzed by considering discrete phase planes, corresponding to the discrete number of open BK channels, and stochastic events correspond to jumps between these planes. The evolution within each plane can be understood from nullclines andlimitcycles,andtheoveralldynamics,e.g.,whether


Introduction
Biological systems are often influenced by discrete, stochastic events, and are therefore appropriately described by hybrid stochastic-deterministic models.Whereas smooth dynamical systems are routinely studied with geometrical techniques, this is still not the case for random dynamical systems such as the ones emerging from hybrid models.The main purpose of the present work is to propose a geometrical method for analyzing hybrid systems, which we illustrate with an application to a hybrid model of cellular electrophysiology.
Many cell types rely on electrical activity to transduce stimuli to signals to be communicated to other cells.In particular, most endocrine cells release hormones as a result of calcium influx via voltage-sensitive Ca 2+ -channels that open during electrical activity, which triggers calcium-dependent exocytosis of hormone-containing vesicles [3,6,22,24,27].Pituitary endocrine cells, such as the prolactin-secreting lactotrophs studied here, are generally small, and hence stochastic ion channel dynamics can have a large influence on the electrical patterns that they exhibit, and hence on the amount of hormone being released [10,25].
Large-conductance potassium (BK) channels play a particular important role in lactotrophs, since they can promote so-called bursting electrical activity, where small action potentials fire from a depolarized plateau [30,32].This effect is biologically relevant, since bursting leads to higher rates of secretion than simpler action potential firing [28,31,32].Previous mathematical deterministic modeling and slow-fast analyses have provided insight into how BK channels promote such bursting [30,33].
However, very few BK channels are active, and mathematical models must therefore take this stochastic and discrete aspect into account.The discrete aspect is particular for BK channels due to their scarcity, and indeed analyses of the effects of stochastic dynamics of other types of ion channels have typically assumed that the system could be described by stochastic differential equations [8, 11-13, 17, 20, 21, 23].
A recent simulation study [25] provided some insight into the role of discrete stochastic BK dynamics for shaping electrical activity in pituitary cells.The authors found numerically that stochastic opening of a BK channel increases the probability of observing a burst if the event happens during the action potential (AP) upstroke or near the AP peak, but lower the burst probability if it occurs after the AP peak.The opposite was observed for stochastic BK channel closing events.
Here we extend the work by Richards et al. [25] by considering a biologically more correct formulation of the control of BK channels [19], and we provide a detailed mathematical analysis of how stochastic opening and closing of discrete BK channels may or may not lead to a burst.Our analysis is based on casting the model in a random-dynamical-system formulation [1] combined with geometric analysis.We note that since the BK channel transition probabilities depend on the membrane potential V , the discrete-noise model that we study is not a so called blinking system where the stochastic dynamics is independent of the deterministic part, i.e., a system purely driven by a discrete random process, e.g., a Markov Chain [2,14].
Our analysis also differs from traditional geometrical analyses of discrete stochastic models of cellular electrical activity, which typically analyze the deterministic model with tools from the theory of (smooth) dynamical systems, and then consider noise as perturbations that "pushes" the system around in phase space [10,25].Here we will show that it is important to understand when and where the "pushes", i.e. the stochastic events, occur.This insight can be obtained by, first, taking advantage of the slow-fast structure of the model, and, then, considering a family of nullclines lying in discrete phase planes corresponding to the discrete number of open BK channels.Together these planes form the fast-subsystem phase space, and stochastic events correspond to jumps between these planes.The dynamics of the system within each plane is determined by geometrical structures such as the nullclines until the next stochastic event, and insight into the overall dynamics can be understood by considering where stochastic jumps between planes occur with respect to the nullclines.

Methods
We devise and analyze a hybrid version of the model by Tabak et al. [30], where overdots indicate differentiation with respect to time t.V is the cellular membrane potential, n is a gating variable for the voltage sensitive K + current (I Kv ), and Ca c is the free cytosolic Ca 2+ concentration.I Ca and I L are (deterministic) voltage dependent Ca 2+ and leak currents, respectively, I SK is the (deterministic) Ca 2+ -gated small-conductance K + (SK) current, while I BK represents the stochastic BK current, which is a function not only of V , but also of the local Ca 2+ concentration, Ca loc , at each BK channel, which depends on the surrounding open Ca 2+ channels (as explained below).C m is the membrane capacitance, α changes current to flux, k c is the Ca 2+ removal rate and f is the ratio of free-to-bound Ca 2+ .The deterministic currents are modeled as where g X represents the whole-cell conductance of the channel X and V X the corresponding reversal potential.The steady-state activation functions, m CaV,∞ and n ∞ , are described by Boltzmann functions, The stochastic current I BK is modeled by where ḡBK is the BK single-channel conductance, n BK the number of BK channels and and for each CaV surrounding the BK channel, where the elements correspond to the transition probabilities between the indicated states in the time interval [t, t + ∆t], provided that ∆t is small.Hence, α and β represent the voltage-dependent Ca 2+ channel opening and closing rates, respectively, and are given by [19,26] In equation (12), k + and k − are the voltage and Ca 2+dependent opening and closing rates for BK channels, and are modeled as in [19], with Ca loc determined by the number of surrounding open CaV channels, n CaV BK,o , with Here, i Ca = ḡCa (V − V Ca ) is the single-channel Ca 2+ current, r the distance between CaVs and the BK channel in a single BK-CaV complex, D Ca the Ca 2+ diffusion constant, [B total ] the total amount of Ca 2+ buffer, and k + B the on-rate of the buffer.When all the surrounding CaVs are closed (i.e., n CaV BK,o = 0), then Ca loc = Ca c .Note that we assume the linear buffer approximation for computing the Ca 2+ profile from n open channels by superimposing n nanodomains found for single, isolated CaVs.
The hybrid system was solved with a fixed time-step procedure implemented in MATLAB/SIMULINK with ∆t = 0.01 ms.A third-order Bogacki-Shampine scheme was used to solve the ODEs (1)- (3).The stochastic part of the model, computing the state of the BK-CaV complex, was updated as follows.At any time point t, a random number ξ uniformly distributed on the interval [0, 1] was generated for each of the surrounding CaV channels, and a transition was made based upon the subinterval in which ξ fell; for example, if the CaV channel was open (O t CaV ) (see the second column of Q CaV defined by ( 13)), it remained open if ξ < 1 − β∆t, otherwise a transition to the closed (C t+∆t CaV ) state occured.Similarly, a random number η uniformly distributed on the interval [0, 1] for the BK channel was generated, and a transition was made based upon the subinterval that η belonged to.
Table 1 reports the parameter values of the model.Matlab code and Simulink schemes implemented for the devised hybrid model with different configurations of the single BK-CaV complex (i.e.different stoichiometries) and their number (i.e.n BK ) are provided as Supplementary Material (see Section "Data availability").
The hybrid model produces different kinds of behavior depending on the configuration of the BK-CaV complexes and their number.Figure 1 shows simulated traces for n BK = 5 BK channels in complexes with 1, 2 or 4 CaVs, which are located either 13 nm or 30 nm from the BK channel of the complex.The membrane voltage V exhibits both single action potential firing as well as so called bursts where small-amplitude action potentials and voltage fluctuations appear from a depolarized plateau.Interesting, bursting appears to be less frequent in the BK-CaV configurations with either many CaVs located close to the BK channel (1:4 stochiometry, r = 13 nm, Fig. 1e), or with few CaVs at a greater distance from the BK channel (1:1 stochiometry, r = 30 nm, Fig. 1b).These two cases correspond respectively to the configurations where we would expect BK channels to open readily or less frequently, as confirmed by the lower traces showing the number of open BK channels.
With more BK channels, bursting seems to be more frequent (Fig. 2, n BK = 15).However, as seen most clearly for the configuration with 1:4 stoichiometry and r = 13 nm (Fig. 2e), the larger number of BK channels tend to make the interburst interval more depolarized (∼ −50 mV) and the active phase of the burst more hyperpolarized (∼ −30 mV), compared to the behavior seen for n BK = 5 (Fig. 1).
To understand the differences seen in the stochastic simulations for the various configurations, we first aim at obtaining a geometric understanding of why the model sometimes produce a burst and sometimes produce a spike.We will then take advantage of this in-sight to explain the behavior seen in the different cases in Figs. 1 and 2.
The time-scale of the gating variable n of the Kvcurrent is τ n = 30 ms, m BK has time scale of a few ms [7,19], whereas Ca c has time scale 1/(f k c ) = 833 ms, which allow us to treat the system as a slow-fast system with Ca c as a slow variable, and the remaining variables (V, n, m BK ) as a hybrid (random) fast subsystem of the full model.
The phase space of this subsystem is composed of 1 + n BK discrete planes, more specifically -since n ∈ [0, 1] -it equals R × [0, 1] × {0, . . ., n BK }.It is useful to consider the V and n nullclines for fixed m BK , the number of open BK channels.The n nullcline is independent of m BK and given by n = n ∞ (V ).The V nullcline, in contrast, depends on m BK and on the (fixed) value of Ca c .
Drawing these nullclines for n BK = 5 and Ca c = 0.4 µM in each of the planes R × [0, 1] × {m BK } shows that for m BK = 0 only a single stable equilibrium exists at V ≈ −15 mV, which is surrounded by an unstable limit cycle and a larger stable limit cycle (Figs.3bc  and 4bc).Increasing m BK , the V nullcline moves downwards and the unstable limit cycle disappears so that the upper stable equilibrium becomes unstable for m BK = 1 and m BK = 2, where three equilibriums are present: the lower one is stable and the one in the middle is a saddle point.For m BK ≥ 3 only the lower stable equilibrium is present.
At the beginning of an action potential when the cell is hyperpolarized, the BK channels tend to close so m BK = 0 most of the time.For this value of m BK the system is attracted to the large limit cycle.However, as V increases the opening rate of the channels increases compared to their closing rate and the fast subsystem   jumps to the planes with m BK > 0, eventually reaching the plane with m BK = 5 (Fig. 3abc).In parallel to the increase in m BK , the Kv gating variable n increases as well, with the result that the fast subsystem is above the V nullcline in the m BK = 5 plane when the trajectory reaches this plane.Consequently, V starts to decrease whereas n still increases.As the cell hyperpolarizes, the BK channels begin to close, eventually followed by a decrease in the n variable.What distinguishes a spike from a burst is the point at which the trajectory reaches the m BK = 0 plane.If it happens above and to the left of the middle part of V nullcline, V will continue to decline, terminating the action potential so that a single spike occurred (Fig. 3).This is true even if a BK channel should open as this would move the V nullcline downwards and the trajectory would remain above.
If, on the other hand, the trajectory reaches the m BK = 0 plane below/to the right of the V nullcline as in Fig. 4c, V starts to increase leading to a second action potential.This scenario might repeat itself several times with the result that a burst is formed.During the burst, the Ca 2+ level Ca c increases (Fig. 4a), which moves the V nullcline downwards (Fig. 4c).This shift increases the probability that the system hits above/to the left of the V nullcline in the m BK = 0 plane, which would terminate the burst.That is, slow feedback from Ca 2+ contributes to controlling the end of the burst.
With less CaVs in the BK-CaVs complexes, or with a greater distance between the CaVs and the BK channel, the steady-state average fraction of open BK channels decreases for any given V (Fig. 5), as expected from the biophysical fact that BK channels are Ca 2+ activated, and hence, if a lower Ca 2+ concentration is present at the BK channel because of fewer or more distant CaVs in its complex, the open probability decreases.This observation implies for example that with only a single CaV located in the BK-CaV complexes, less BK channels will open as the cell depolarizes during the beginning of the action potential, compared to the case with 1:4 stoichiometry.For example, in Fig. 6 (1:1 stoichiometry and r = 13 nm), only 1-2 BK channels are open during most of the first increase in V (red points), in contrast to the 2-4 BK channels being open during the upstroke in Fig. 4 (1:4 stoichiometry and r = 13 nm).This implies that the system is below the V nullcline for a longer time, leading to a first peak at V ∼ 0 mV, compared to the peak at V ∼ −10 mV for the case of 1:4 stoichiometry.Moreover, during the beginning of the downstroke (blue points in Fig. 6c) the system is relatively close to the V nullcline, meaning that V does not decrease rapidly.Altogether, the long time that the systems spends at very depolarized V , allows the variable n to increase more than in the previous case.The result is that when BK channels eventually close and the system reaches the m BK = 0 plane, it will often fall within the unstable limit cycle and begin to spiral counter-clockwise.Even if one or two BK channels open, the system is still in the region of the (V, n) plane where the unstable limit cycle is lying in the m BK = 0 plane.Only when Ca c has increased sufficiently, so that the V nullcline has moved sufficiently downwards, is the system able to escape and terminate the burst.To obtain spiking, the system must fall outside and to the left of the unstable limit cycle when reaching the m BK = 0 plane, which is confirmed by our simulations (not shown).
To summarize, compared to the 1:4 stoichiometry scenario, the fact that the opening probability for BK channels is lower causes n to increase so much that the system often falls on the inside of the unstable limit cycle when returning to the m BK = 0 plane.This mechanism explains why bursting is more frequent with fewer CaVs in the BK-CaV complexes for r = 13 nm, compare panels (a) and (e) in Fig. 1.
A similar explanation underlies the increased propensity for bursting when CaVs are located 30 nm (rather than 13 nm) from the BK channel in the complex with 1:4 stoichiometry (panels (e) and (f) in Fig. 1).Again, BK channels tend to open less during the upstroke so that n can increase more, which eventually makes V decrease.As V decreases, the BK channels tend to close and the system reaches the m BK = 0 plane at fairly high n values (Fig. 7) so that the system might fall inside or very close to the unstable limit cycle, but to the right of the V nullcline, thus creating small-amplitude oscillations (Fig. 7).This mechanism is similar to the one described producing bursting for r = 13 nm and 1:4 stoichiometry, but it is more frequent for r = 30 nm since n typically will be higher at the action potential peaks, compare the red traces in panels (e) and (f) in Fig. 1.
With more BK channels in the cells, a higher number of BK channels will be open in spite of the same opening probability.For example, with n BK = 15 instead of nBK = 5, three times as many BK channels will be open on average, in spite of all other configurations (membrane potential V , number of CaVs per BK-CaV complex, and CaV-to-BK distance r) being identical.Thus, during the upstroke, the system quickly reaches m BK > 5, and the many open BK channels stop the increase in V and lead to a first peak at V < −20 mV (Fig. 8).As a consequence of the relatively weak depolarization of the membrane potential, few Kv channels become activated, i.e., n remains low (n < 0.08 For m BK = 0 the trajectories starting outside the unstable limit cycle converge to a stable period orbit (gray).The simulation in panel (a) is projected onto the plane with the corresponding value of m BK using dots with colors as in panel (a).Note how the system returns to the m BK = 0 plane to the left of the V nullcline (green points), which forces V to decrease further, ending the action potential.(c) Fig. 4 An example of a burst for n BK = 5 with 1:4 stoichiometry and r = 13 nm.It corresponds to the burst of Fig. 1e for 4.12 ≤ t ≤ 4.35.Legends as in Fig. 3.In the first subplot of panel (c), for m BK = 0, the V -nullcline (light-blue) and the unstable limit cycle (violet) are calculated for Cac = 0.4 µM (solid), Cac = 0.42 µM (dashed) and Cac = 0.46 µM (dotted; the limit cycle is very small with center (-17 mV, 0.23)).Note how the system returns to the m BK = 0 plane to the right of the V -nullclines the first two times (first time, green points to the right of the solid light-blue line; second time, black points to the right of dashed light-blue line), leading to a new increase of V .Only the third time, when Cac has increased, moving the V nullcline further downwards (dotted light-blue line), does the system (brown points) return to the m BK = 0 plane to the left of the V nullcline, which causes a further decrease of V , terminating the burst.
in Fig. 8).Geometrically, this means that the system remains in the lower part of the (V, n) planes.Consequently, when the BK channels close during the downstroke and the system reaches the m BK = 0 plane, it will be to the right of the V nullcline (Fig. 8c), and the system will produce additional oscillations going from V ≈ −40 mV to V ≈ −25 mV, until Ca c increases sufficiently to move the V nullcline downward so that the system falls in the region above the V nullcline.
In contrast to the previous cases, this occurs for much lower n values, and hence the system does not follow the trajectories that reach a minimum for V at −60 mV.

Discussion
Understanding how complex dynamics arises in nonlinear, hybrid deterministic-stochastic models is important for providing insight into the role of underlying mechanisms and for controlling the corresponding systems.Geometrical methods have proven highly useful for deterministic systems, also in biology [9,15,16].Whereas deterministic models driven by continuous stochastic processes, e.g., Wiener processes, have been quite extensively studied [4,8,20], this is not the case for discrete hybrid systems.However, some results do exist for systems driven by a Markov chain independent of the other (deterministic) variables, so called blinking systems [2,14].The model of electrical activity in a pituitary cell with stochastic ion channels that we studied here does not fall within either of the examples mentioned above, but is a truly random dynamical system [1], since the transition rates of the stochastic variable (m BK ) depend on the deterministic variables (V ).Recent studies have studied similar hybrid models of stochastic electrical activity by analyzing the corresponding "average" deterministic model with geometric tools, and then interpreting ion channel noise as random perturbations ("pushes") of the deterministic system [10,25].
In contrast, here we have shown how one can work directly with the hybrid system.This is achieved by analyzing the hybrid (fast subsystem) phase space as the union of discrete planes, each one corresponding to a certain value of the discrete stochastic variable m BK indicating the number of open BK channels.In the terminology of random dynamical systems, each of these planes is a fiber, and the random part of the model corresponds to jumps between these fibers.We showed that the locations of geometric structures, in particular nullclines and an unstable limit cycle, govern the behavior for fixed m BK , and that the overall dynamics, e.g., whether the model produces a spike or a burst, is determined by the location at which the system jumps from one plane to another, in particular, the point at which the system reaches the m BK = 0 plane plays an important role.
To reach this description, we took advantage of the slow-fast structure of the model.Since the Ca 2+ variable Ca c operates on a slower timescale than the other variables, it can be treated as a (slowly varying) parameter in the fast subsystem.For our model, this assumption has the big advantage that the fibers corresponding to fixed m BK values become two-dimensional, which helps geometrical reasoning.The slow dynamics of Ca c was taken into account by considering how the relevant geometrical structures, in particular the V -nullcline, move as Ca c changes.
The strength of our approach is maybe best exemplified by its ability to explain counter-intuitive numerical results.We noted that spiking is more frequently seen both when BK channels are located in complexes with four CaV channels at a distance of r = 13 nm (Fig. 1e), which would lead to a high BK open probability, and when BK channels are associated with just a single CaV at a large distance (r = 30 nm, Fig. 1c), corresponding to a low BK open probability.Thus, there seems to be a window of intermediate BK opening probabilities where bursting is favored.We explained this by noticing that if the BK channels open to readily during the upstroke, the gating variable n does not increase very much since the V variable begins to decrease early.When returning to the m BK = 0 plane the system is often below that unstable limit cycle, but at n large enough to be above the V nullclines, which ends the action potential.If the BK channels open more slowly, V and consequently n increase more, so that the system returns to the m BK = 0 plane near or even inside the unstable limit cycle, thus creating a second action potential, i.e., a burst, before hyperpolarizing.In contrast, if the BK channel do not open sufficiently, so that m BK ≤ 1 except on rare occations (see Fig. 1c), the system basically follows the stable limit cycle, corresponding to spiking, in the m BK = 0 plane, and the large orbit in the m BK = 1 plane (see gray curves in, e.g., Fig. 3c).Biologically, this suggests that BK-CaV complexes could be appropriately tuned to increase the burst frequency.
Our analysis can also explain geometrically the observations by Richards et al. [25] that opening a BK channel just before the action potential peak increases (c) Fig. 6 An example of a burst for n BK = 5 with 1:1 stoichiometry and r = 13 nm.It corresponds to the burst of Fig. 1a for 4.82 ≤ t ≤ 4.96.Legends as in Fig. 3.Note how the system returns to the m BK = 0 plane (first subplot of panel (c)) within the unstable limit cycle (green points inside the violet oval) leading to a new increase of V (magenta points).Only the second time, when Cac has increased moving the V nullcline slightly downwards, does the system return to the m BK = 0 plane to the left of the V nullcline (cyan points to the left of the dotted light-blue line with Cac = 0.44 µM), which causes a further decrease of V , terminating the burst.For Cac = 0.44 µM, the unstable limit cycle is also reduced.
the probability of observing a burst, whereas opening a BK channel during the downstroke of the action potential reduces the chance of producing a burst.If a BK channel opens at high V , the V nullcline moves down and V will stop increasing and begin decreasing, which leads to less activation of Kv channels (smaller n).The results is that the system returns to the m BK = 0 plane at lower n values, where it is more likely to fall below and on the right of the middle branch of the V nullcline.
In contrast, if a BK channel opens during the down-stroke when both V and n are decreasing, the downward shift of the V nullcline will accelerate the decrease of V , so that when the system eventually reaches the m BK = 0 plane, it will more likely hit above and to the left of the middle part of the V nullcline.
In summary, we have presented a -to the best of our knowledge -novel method that combines ideas from standard geometrical analysis of smooth dynamical systems with a picture taken from random dynamical systems, where stochastic events correspond to random jumps between fibers.This approach, which successfully allowed us to explain complex behavior in a hybrid model of electrical activity influenced by stochastic ion channel dynamics in a pituitary cell, should be useful for similar models in cell biology as well as for other applications of nonlinear, hybrid models.

Statements and Declarations
Funding: The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Fig. 1
Fig.1The cellular membrane potential V (upper plot in each panel), the gating variable n for the Kv current (second plot), the free cytosolic Ca 2+ concentration Cac (third plot) and the number of open BK channels (lower plot, m BK open) with respect to time by assuming a number of total BK (n BK ) equal to 5. For each BKCa-CaV ion channel complex, different stoichiometries are employed: 1:1 for the first row (panels (a) and (b)); 1:2 for the second row (panels (c) and (d)); 1:4 for the third row (panels (e) and (f)).Also, two different values for the distance r between CaVs and BK channel in a single BK-CaV complex are considered: r = 13 nm for the first column (panels (a), (c) and (e)); r = 30 nm for the second column (panels (b), (d) and (f)).

Fig. 2 V
Fig. 2 V, n, Cac and m BK as functions of time for n BK = 15.Details as in Fig. 1.

Fig. 3
Fig. 3 Single action potential example.A single action potential (spike) for n BK = 5 BK channels in complexes with 4 CaVs (i.e., 1:4 stoichiometry) located 13 nm (= r) from the BK of the complex.(a) V , n, Cac and m BK as functions of time.It corresponds to the spike of Fig. 1e for 4.64 ≤ t ≤ 4.72.The colors of the curves indicate different phases of the action potential for easier comparison to panels (b) and (c).(b) Projection of the simulation in panel (a) onto the phase space of the fast (V, n, m BK ) subsystem (colors as in panel (a)).The phase space is composed of the gray planes given by constant m BK , and random opening and closing of BK channels correspond to jumps between these planes.For fixed m BK , the V (light blue) and n (red) nullclines are shown for Cac = 0.4 µM.(c) Phase planes for fixed m BK as indicated and Cac = 0.4 µM, corresponding to the gray planes in panel (b), with V (light blue) and n (red) nullclines.The red circle for m BK = 0 represent the first sample of the time simulation.The violet oval for m BK = 0 indicates the unstable limit cycle.Gray asterisks indicate initial conditions for different deterministic orbits obtained by keeping m BK fixed and equal to the value of the corresponding panel.For m BK = 0 the trajectories starting outside the unstable limit cycle converge to a stable period orbit (gray).The simulation in panel (a) is projected onto the plane with the corresponding value of m BK using dots with colors as in panel (a).Note how the system returns to the m BK = 0 plane to the left of the V nullcline (green points), which forces V to decrease further, ending the action potential.

Fig. 7
Fig.7An example of a burst for n BK = 5 with 1:4 stoichiometry and r = 30 nm.It corresponds to the burst of Fig.1ffor 3.68 ≤ t ≤ 3.91.Legends as in Fig.3.Note how the system returns to the m BK = 0 plane (first subplot of panel (c)) to the right of the V -nullclines the first two times (first time, green points to the right of the solid light blue line for Cac = 0.4; second time, black points to the right of dashed blue line with Cac = 0.45 µM), leading to a new increase of V .Only the third time, when Cac has increased, moving the V nullcline downwards, does the system return to the m BK = 0 plane to the left of the V nullcline (brown points to the left of the dotted blue line with Cac = 0.49 µM), which causes a further decrease of V terminating the burst.Note how the oval unstable limit cycle for m BK = 0 reduces by increasing Cac, approaching to the infinitesimal cycles with centres (-17 mV, 0.24) and (-18 mV, 0.23) for Cac = 0.45 and Cac = 0.49 µM, respectively.

Table 1
Parameter values of the pituitary model.