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

Abstract

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 propose 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 and limit cycles, and 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. Our approach is generally applicable to other scenarios to study discrete random dynamical systems defined by hybrid stochastic–deterministic models.

1 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²⁺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, therefore, 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,12,13, 17, 20, 21, 23]. Such a continuous-noise approach can be justified through an appropriate approximation to the underlying Markov chain dynamics of the ion channel population in the case of a sufficiently high number of active ion channels [11,12,13].

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 phase space of the fast subsystem, 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.

2 Methods

We devise and analyze a hybrid version of the model by Tabak et al. [30],

$$\begin{aligned} \dot{V}&= -\big (I_\textrm{Ca}(V)+I_\textrm{Kv}(V,n)+I_\textrm{SK}(V,{\textrm{Ca}_c}) \nonumber \\&\quad +I_\textrm{BK}(V,{{\textrm{Ca}_\textrm{loc}}})+I_{L}(V)\big )/C_m, \end{aligned}$$

(1)

$$\begin{aligned} \dot{n}&= (n_{\infty }(V)-n)/{\tau _n}, \end{aligned}$$

(2)

$$\begin{aligned} \dot{\textrm{Ca}_c}&= f (\alpha \, I_\textrm{Ca}(V) - k_c\, {\textrm{Ca}_c}), \end{aligned}$$

(3)

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_\textrm{Kv}$), and ${\textrm{Ca}_c}$ is the free cytosolic Ca²⁺ concentration. $I_\textrm{Ca}$ and $I_{L}$ are (deterministic) voltage-dependent Ca²⁺ and leak currents, respectively, $I_\textrm{SK}$ is the (deterministic) Ca²⁺-gated small-conductance K⁺ (SK) current, while $I_\textrm{BK}$ represents the stochastic BK current, which is a function not only of V, but also of the local Ca²⁺ concentration, ${\textrm{Ca}_\textrm{loc}}$, at each BK channel, which depends on the surrounding open Ca²⁺ channels (as explained below). $C_m$ is the membrane capacitance, $\alpha$ changes current to flux, $k_c$ is the Ca²⁺ removal rate and f is the ratio of free-to-bound Ca²⁺.

The deterministic currents are modeled as

$$\begin{aligned} I_\textrm{Ca}(V)&=g_\textrm{Ca}\,m_{\textrm{CaV},\infty }(V)(V-V_\textrm{Ca}), \end{aligned}$$

(4)

$$\begin{aligned} I_\textrm{Kv}(V,n)&=g_K\,n\,(V-V_K), \end{aligned}$$

(5)

$$\begin{aligned} I_\textrm{SK}(V,\textrm{Ca}_c)&=g_\textrm{SK}\,s_{\infty }(\textrm{Ca}_c)(V-V_{K}), \end{aligned}$$

(6)

$$\begin{aligned} I_{L}(V)&=g_{L}\,(V-V_{L}), \end{aligned}$$

(7)

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_{\textrm{CaV},\infty }$ and $n_{\infty }$, are described by Boltzmann functions,

$$\begin{aligned} m_{\textrm{CaV},\infty }(V)&=\frac{1}{1+ \exp ((v_{m}-V)/s_m)}, \end{aligned}$$

(8)

$$\begin{aligned} n_{\infty }(V)&=\frac{1}{1+ \exp ((v_{n}-V)/s_n)}, \end{aligned}$$

(9)

and the steady-state calcium-dependent activation function, $s_{\infty }$, by

$$\begin{aligned} s_{\infty }(\textrm{Ca}_c)=\frac{\textrm{Ca}_c^{2}}{\textrm{Ca}_c^2+k_s^2}. \end{aligned}$$

(10)

The randomness of the system is due to the stochastic current $I_\textrm{BK}$, which is modeled by

$$\begin{aligned} \begin{aligned} I_\textrm{BK}&= \sum _{i=1}^{n_\textrm{BK}} {{\bar{g}}_\textrm{BK} {\textbf{1}}_{O_{\textrm{BK},i}}} (V-V_K) \\&={{\bar{g}}_\textrm{BK} m_\textrm{BK}(V-V_K)}, \end{aligned} \end{aligned}$$

(11)

where ${{\bar{g}}_\textrm{BK}}$ is the BK single-channel conductance, $n_\textrm{BK}$ the number of BK channels and ${\textbf{1}}_{O_{\textrm{BK},i}}$ an index function that equals 1 if the i-th BK is open and 0 otherwise, so that $m_\textrm{BK}$ indicates the total number of open BK channels. In order to compute the state of the BK channels, we exploit a model of single-channel gating with two states (closed and open), whose dynamics are regulated globally via membrane potential V, and locally via the Ca²⁺ nanodomains below the mouth of stochastic CaV channels surrounding the single BK channel [18, 19], forming an ion channel BK-CaV complex with a stoichiometry of 1–4 CaV channels per BK channel [5, 29]. Therefore, in order to model the stochastic gating of BK channels, we describe the transition from one state (closed or open) to another for each BK channel and the surrounding CaVs. ($C_X^t$ corresponds to the closed state and $O_X^t$ to the open state of the channel X at time t.) Then, we define the transition matrix for the single BK,

$$\begin{aligned} Q_\textrm{BK}&=\begin{bmatrix} P\left( C_\textrm{BK}^{t+ \varDelta t} |C_\textrm{BK}^t\right) &{}\quad P\left( C_\textrm{BK}^{ t+ \varDelta t} |O_\textrm{BK}^t\right) \\ P\left( O_\textrm{BK}^{t+ \varDelta t} |C_\textrm{BK}^t\right) &{}\quad P\left( O_\textrm{BK}^{t+ \varDelta t} |O_\textrm{BK}^t\right) \end{bmatrix}\nonumber \\&=\begin{bmatrix} 1-k^+ \varDelta t &{}\quad k^- \varDelta t \\ k^+ \varDelta t&{}\quad 1-k^- \varDelta t \\ \end{bmatrix}, \end{aligned}$$

(12)

and for each CaV surrounding the BK channel,

$$\begin{aligned} Q_\textrm{CaV}&= \begin{bmatrix} P\left( C_\textrm{CaV}^{t+ \varDelta t} |C_\textrm{CaV}^t\right) &{}\quad P\left( C_\textrm{CaV}^{t+ \varDelta t} |O_\textrm{CaV}^t\right) \\ P\left( O_\textrm{CaV}^{t+ \varDelta t} |C_\textrm{CaV}^t\right) &{}\quad P\left( O_\textrm{CaV}^{t+ \varDelta t} |O_\textrm{CaV}^t\right) \end{bmatrix}\nonumber \\&=\begin{bmatrix} 1-\alpha \varDelta t &{}\quad \beta \varDelta t \\ \alpha \varDelta t&{}\quad 1-\beta \varDelta t \end{bmatrix}, \end{aligned}$$

(13)

where the elements correspond to the transition probabilities between the indicated states in the time interval $[t, t+\varDelta t]$, provided that $\varDelta t$ is small. Hence, $\alpha$ and $\beta$ represent the voltage-dependent Ca²⁺ channel opening and closing rates, respectively, and are given by [19, 26]

$$\begin{aligned} \alpha (V)&= \frac{m_{\textrm{CaV},\infty }(V)}{\tau _\textrm{CaV}}, \end{aligned}$$

(14)

$$\begin{aligned} \beta (V)&= \frac{1-m_{\textrm{CaV},\infty }(V)}{\tau _\textrm{CaV}}. \end{aligned}$$

(15)

In equation (12), $k^+$ and $k^-$ are the voltage and Ca²⁺-dependent opening and closing rates for BK channels, and are modeled as in [19],

$$\begin{aligned} k^+(V,\textrm{Ca}_\textrm{loc})&= w^+(V)f^+(\textrm{Ca}_\textrm{loc})\nonumber \\&= w_0^+e^{-w_\textrm{co}V} \frac{1}{1+\left( \frac{K_\textrm{co}}{{\textrm{Ca}_\textrm{loc}}}\right) ^{n_\textrm{co}}}, \end{aligned}$$

(16)

$$\begin{aligned} k^-(V,\textrm{Ca}_\textrm{loc})&= w^-(V)f^-(\textrm{Ca}_\textrm{loc}) \nonumber \\&= w_0^-e^{-w_\textrm{oc}V} \frac{1}{1+\left( \frac{\textrm{Ca}_\textrm{loc}}{{K_\textrm{oc}}}\right) ^{n_\textrm{oc}}}, \end{aligned}$$

(17)

Table 1 Parameter values of the pituitary model

with Ca$_\textrm{loc}$ determined by the number of surrounding open CaV channels, $n_{\textrm{CaV}_{\textrm{BK},o}}$,

$$\begin{aligned} \textrm{Ca}_\textrm{loc}={n_{\textrm{CaV}_{\textrm{BK},o}}} \textrm{Ca}_{o} +\textrm{Ca}_c, \end{aligned}$$

(18)

with

$$\begin{aligned} \textrm{Ca}_o = \frac{i_\textrm{Ca}}{8 \pi r D_\textrm{Ca}F}\exp {\left[ \frac{-r}{\sqrt{\frac{D_\textrm{Ca}}{k^+_B[B_\textrm{total}]}}}\right] }. \end{aligned}$$

(19)

Here, $i_\textrm{Ca} = {\bar{g}}_\textrm{Ca}(V-V_\textrm{Ca})$ is the single-channel Ca²⁺ current, r the distance between CaVs and the BK channel in a single BK-CaV complex, $D_\textrm{Ca}$ the Ca²⁺ diffusion constant, $[B_\textrm{total}]$ the total amount of Ca²⁺ buffer and $k_B^+$ the rate of the buffer. When all the surrounding CaVs are closed (i.e., $n_{\textrm{CaV}_{\textrm{BK},o}}=0$), then Ca$_\textrm{loc}=\textrm{Ca}_c$. Note that we assume the linear buffer approximation for computing the Ca²⁺ 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 $\varDelta t =0.01$ ms. Since we are interested in long-term behavior, the choice of the initial values is irrelevant. 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 $\xi$ 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 $\xi$ fell; for example, if the CaV channel was open ($O^t_\textrm{CaV}$) (see the second column of $Q_\textrm{CaV}$ defined by (13)), it remained open if $\xi <1-\beta \varDelta t$; otherwise, a transition to the closed ($C^{t+\varDelta t}_\textrm{CaV}$) state occurred. Similarly, a random number $\eta$ uniformly distributed on the interval [0, 1] for the BK channel was generated, and a transition was made based upon the subinterval that $\eta$ 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_\textrm{BK}$) are provided in a freely accessible online repository (see Section “Data availability”).

Fig. 1

Cellular membrane potential V (upper plot in each panel), the gating variable n for the $K_v$ current (second plot), the free cytosolic Ca²⁺ concentration Ca$_c$ (third plot) and the number of open BK channels (lower plot, $m_\textrm{BK}$ open) with respect to time by assuming a number of total BK ($n_\textrm{BK}$) equal to 5. For each BKCa-CaV ion channel complex, different stoichiometries are employed: 1:1 for the first row (a and b); 1:2 for the second row (c and d); 1:4 for the third row (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 (a, c and e); $r=30$ nm for the second column (b, d and f)

3 Results

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_\textrm{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 stoichiometry, $r=13$ nm, Fig. 1e), or with few CaVs at a greater distance from the BK channel (1:1 stoichiometry, $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_\textrm{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 ($\sim -\,50$ mV) and the active phase of the burst more hyperpolarized ($\sim -\,30$ mV), compared to the behavior seen for $n_\textrm{BK}=5$ (Fig. 1).

Fig. 2

$V,n,\textrm{Ca}_c$ and $m_\textrm{BK}$ as functions of time for $n_\textrm{BK}=15$. Details as in 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 produces a burst and sometimes produces a spike. We will then take advantage of this insight to explain the behavior seen in the different cases in Figs. 1 and 2.

The timescale of the gating variable n of the Kv current is $\tau _n=30$ ms, $m_\textrm{BK}$ has timescale of a few ms [7, 19], whereas Ca$_c$ has timescale $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_\textrm{BK})$ as a hybrid (random) fast subsystem of the full model.

The phase space of this subsystem is composed of $1+n_\textrm{BK}$ discrete planes, more specifically—since $n\in [0,1]$—it equals ${\mathbb {R}}\times [0,1]\times \{0,\ldots ,n_\textrm{BK}\}$. It is useful to consider the V and n nullclines for fixed $m_\textrm{BK}$, the number of open BK channels. The n nullcline is independent of $m_\textrm{BK}$ and given by $n=n_{\infty }(V)$. The V nullcline, in contrast, depends on $m_\textrm{BK}$ and on the (fixed) value of Ca$_c$.

Drawing these nullclines for $n_\textrm{BK}=5$ and Ca$_c=0.4\,\upmu$M in each of the planes ${\mathbb {R}}\times [0,1]\times \{m_\textrm{BK}\}$ shows that for $m_\textrm{BK}=0$ only a single stable equilibrium exists at $V\approx -15$ mV, which is surrounded by an unstable limit cycle and a larger stable limit cycle (Figs. 3b, c and 4b, c). Increasing $m_\textrm{BK}$, the V nullcline moves downwards and the unstable limit cycle disappears so that the upper stable equilibrium becomes unstable for ${m_\textrm{BK}=1}$ and $m_\textrm{BK}=2$, where three equilibriums are present: The lower one is stable and the one in the middle is a saddle point. For $m_\textrm{BK}\ge 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_\textrm{BK}=0$ most of the time. For this value of $m_\textrm{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_\textrm{BK}>0$, eventually reaching the plane with $m_\textrm{BK}=5$ (Fig. 3a–c). In parallel with the increase in $m_\textrm{BK}$, the Kv gating variable n increases as well, with the result that the fast subsystem is above the V nullcline in the $m_\textrm{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_\textrm{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_\textrm{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²⁺ 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_\textrm{BK}=0$ plane, which would terminate the burst. That is, slow feedback from Ca²⁺ 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²⁺ activated, and hence, if a lower Ca²⁺ concentration is present at the BK channel because of fewer or more distant CaVs in its complex, the open probability decreases.

Fig. 3

Single action potential example. A single action potential (spike) for $n_\textrm{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, Ca$_c$ and $m_\textrm{BK}$ as functions of time. It corresponds to the spike of Fig. 1e for $4.64 \le t \le 4.72$. The colors of the curves indicate different phases of the action potential for easier comparison to (b) and (c). b Projection of the simulation in (a) onto the phase space of the fast $(V,n,m_\textrm{BK})$ subsystem (colors as in (a)). The phase space is composed of the gray planes given by constant $m_\textrm{BK}$, and random opening and closing of BK channels correspond to jumps between these planes. For fixed $m_\textrm{BK}$, the V (light blue) and n (red) nullclines are shown for Ca$_c=0.4\,\upmu$M. c Phase planes for fixed $m_\textrm{BK}$ as indicated and Ca$_c=0.4\,\upmu$M, corresponding to the gray planes in (b), with V (light blue) and n (red) nullclines. The red circle for $m_\textrm{BK}=0$ represent the first sample of the time simulation. The violet oval for $m_\textrm{BK}=0$ indicates the unstable limit cycle. Gray asterisks indicate initial conditions for different deterministic orbits obtained by keeping $m_\textrm{BK}$ fixed and equal to the value of the corresponding panel. For $m_\textrm{BK}=0$, the trajectories starting outside the unstable limit cycle converge to a stable period orbit (gray). The simulation in (a) is projected onto the plane with the corresponding value of $m_\textrm{BK}$ using dots with colors as in (a). Note how the system returns to the $m_\textrm{BK}=0$ plane to the left of the V nullcline (green points), which forces V to decrease further, ending the action potential. (Color figure online)

Fig. 4

An example of a burst for $n_\textrm{BK}=5$ with 1:4 stoichiometry and $r=13$ nm. It corresponds to the burst of Fig. 1e for $4.12 \le t \le 4.35$. Legends as in Fig. 3. In the first subplot of (c), for $m_\textrm{BK}=0$, the V-nullcline (light blue) and the unstable limit cycle (violet) are calculated for Ca$_c=0.4\,\upmu$M (solid), Ca$_c=0.42\,\upmu$M (dashed) and Ca$_c=0.46\,\upmu$M (dotted; the limit cycle is very small with center (− 17 mV, 0.23)). Note how the system returns to the $m_\textrm{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 Ca$_c$ has increased, moving the V nullcline further downwards (dotted light blue line), does the system (brown points) return to the $m_\textrm{BK}=0$ plane to the left of the V nullcline, which causes a further decrease of V, terminating the burst. (Color figure online)

Fig. 5

Steady-state BK activation function, $m_{\textrm{BK},\infty }$ (see Eq. 29 in [19])

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\sim 0$ mV, compared to the peak at $V\sim -\,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_\textrm{BK}=0$ plane, it will often fall within the unstable limit cycle and begin to spiral counterclockwise. 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_\textrm{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_\textrm{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_\textrm{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_\textrm{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.

Fig. 6

An example of a burst for $n_\textrm{BK}=5$ with 1:1 stoichiometry and $r=13$ nm. It corresponds to the burst of Fig. 1a for $4.82 \le t \le 4.96$. Legends as in Fig. 3. Note how the system returns to the $m_\textrm{BK}=0$ plane (first subplot of (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 Ca$_c$ has increased moving the V nullcline slightly downwards, does the system return to the $m_\textrm{BK}=0$ plane to the left of the V nullcline (cyan points to the left of the dotted light blue line with Ca$_c=0.44\,\upmu$M), which causes a further decrease of V, terminating the burst. For Ca$_c=0.44\,\upmu$M, the unstable limit cycle is also reduced. (Color figure online)

Fig. 7

An example of a burst for $n_\textrm{BK}=5$ with 1:4 stoichiometry and $r=30$ nm. It corresponds to the burst of Fig. 1f for $3.68 \le t \le 3.91$. Legends as in Fig. 3. Note how the system returns to the $m_\textrm{BK}=0$ plane (first subplot of (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 Ca$_c=0.4$; second time, black points to the right of dashed blue line with Ca$_c=0.45\,\upmu$M), leading to a new increase of V. Only the third time, when Ca$_c$ has increased, moving the V nullcline downwards, does the system return to the $m_\textrm{BK}=0$ plane to the left of the V nullcline (brown points to the left of the dotted blue line with Ca$_c=0.49\,\upmu$M), which causes a further decrease of V terminating the burst. Note how the oval unstable limit cycle for $m_\textrm{BK}=0$ reduces by increasing Ca$_c$, approaching to the infinitesimal cycles with centers (− 17 mV, 0.24) and (− 18 mV, 0.23) for Ca$_c=0.45$ and Ca$_c=0.49\,\upmu$M, respectively. (Color figure online)

Fig. 8

An example of a burst for $n_\textrm{BK}=$ 15 with 1:4 stoichiometry and $r=13$ nm. It corresponds to the burst of Fig. 2e for $0.93 \le t \le 1.16$. Legends as in Fig. 3. Note how the system returns to the $m_\textrm{BK}=0$ plane (first subplot of (c) and the relative zoom in reported in the second subplot) to the right of the V-nullclines more times and within the unstable limit cycle (see magenta, cyan, black, yellow, green, brown and violet points to the right of the solid light blue line for Ca$_c \approx 0.38$), leading to additional oscillations of V from − 40 to − 25 mV. Only, when Ca$_c$ has increased, moving the V nullcline downwards, does the system return to the $m_\textrm{BK}=0$ plane to the left of the V nullcline (orange points to the left of the dotted blue line with Ca$_c=0.41\,\upmu$M), which causes a further decrease of V terminating the burst. Note how the oval unstable limit cycles for $m_\textrm{BK}=0$ reduce by increasing Ca$_c$: compare the greater violet oval for Ca$_c=0.38\,\upmu$M with the smaller one for Ca$_c=0.41\,\upmu$M. (Color figure online)

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_\textrm{BK}=15$ instead of $n\textrm{BK}=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_\textrm{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$ 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_\textrm{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\approx -\,40$ mV to $V\approx -\,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.

4 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 proved 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_\textrm{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_\textrm{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_\textrm{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_\textrm{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²⁺ 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_\textrm{BK}$ values become two-dimensional (the (V, n) phase planes shown in panels (b) and (c) in Figs. 3, 4, 6, 7 and 8), 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 counterintuitive 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_\textrm{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_\textrm{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_\textrm{BK}\le 1$ except on rare occasions (see Fig. 1c), the system basically follows the stable limit cycle, corresponding to spiking, in the $m_\textrm{BK}=0$ plane, and the large orbit in the $m_\textrm{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 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_\textrm{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 downstroke 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_\textrm{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. We conclude that, similar to deterministic systems, the use of geometric slow–fast analysis of stochastic models allows analysis of the dynamics of lower-dimensional subsystems, and the overall model behavior can be understood globally by appropriate integration of the information obtained for the subsystems. Future work should investigate more theoretical aspects of our analytical approach for discrete-noise random dynamical systems, as well as aim to develop the proposed method further.