Review of the deterministic model

We base our pituitary cell model on the model by Tabak et al. [16], a type of Hodgkin-Huxley model [27], where the membrane potential V is controlled by five currents: an L-type Ca²⁺ current (I_Ca), three K⁺ currents (I_K, I_SK and I_BK) and a leak current (I_leak). Then V is described by (1) where C is the membrane capacitance. Note that, since we implement noise by directly considering the opening and closing of individual channels, we have no need for an explicit noise current. The Ca, K, and BK channels are assumed to be voltage-gated, whereas SK is Ca²⁺-gated. Although BK channels are known to be also Ca²⁺-gated, BK channels are usually adjacent to Ca²⁺ channels and so controlled by calcium microdomains [28] meaning that, since the Ca²⁺ concentration at BK channels then attains equilibrium in microseconds, BK channel activation can effectively be modelled as just voltage-gated [29, 30]. The difference between the K and BK channels lies in the single channel conductance (g_1,BK is twenty times larger than g_1,K), the number of channels per cell (with hundreds of K channels and only ten or fewer BK channels), and the time constant (τ_n is six times longer than τ_BK).

Activation of the Ca, K, SK, and BK channels are described by gating variables m, n, s and f respectively, each of which should be thought of as the probability that a given channel is open. These types of model often reduce the number of variables by replacing one or more gating variables by steady-state activation functions. However, since we wish to model realistic noise in all channels, we cannot do this. The Ca²⁺ and K⁺ currents are then given by (2) and the leak current by I_leak = g_l(V − V_l). Here V_Y is the reversal potential for ion Y and g_X the maximal conductance for channel X. Here and in the following we will use to X to refer to an unspecified channel type (either Ca, K, SK or BK). The conductances are related to the single channel conductances via g_X = g_1,X N_X, where N_X is the number of a given channel type.

The dynamics of the gating variables are governed by (3) where τ_X are the time constants (assumed to be voltage-independent). The steady-state activation functions are given by (4) with k_s, v_X and s_X the calcium-midpoint, voltage-midpoint and slope parameters respectively. The dependence of s_∞ on [Ca] is as in previous models of both endocrine and pancreatic β-cells [18, 31].

The final model variable, the intracellular free Ca²⁺ ion concentration [Ca], evolves according to (5) where the first term describes flux through Ca²⁺ channels and the second term captures calcium extrusion. Here α is the conversion factor from charge to concentration, k_c is the extrusion rate, and f_c is the fraction of intracellular calcium present as free Ca²⁺ ions. Note that f_c should not be confused with the BK gating variable f. In total, this leads to a six variable model: the membrane potential (V), the four gating variables (m, n, s, f), and the calcium concentration ([Ca]). For parameter values, initial conditions and how we solve this system numerically, see the Methods and S1 Text.

The stochastic model: Adding ion channel noise

To understand the role of channel noise in pituitary cells, we determine the gating variables via a stochastic process. For each channel type this involves keeping track of the number of open channels. For concreteness, for the remainder of this section, consider just BK channels (described by f); the treatment of the other channel types is similar. The equation for f in Eq (3) can be rewritten (6) where β(V) = f_∞(V)/τ_BK and γ(V) = (1 − f_∞(V))/τ_BK. In this form, it is clear that this can be reinterpreted as a master equation with f the probability P_op that a given single BK channel is open and (1 − f) the probability P_cl it is closed: (7) Seen as a two-state Markov process (with non-constant transition rates), we then have, for each channel, (8)

At any given time, N_op channels are open and N_cl are closed, with N_BK = N_op + N_cl the total (fixed) number of channels and the gating variable f given by f = N_op/N_BK. In a short time step Δt, β(V)Δt is approximately the probability that a single closed channel opens and γ(V)Δt the probability that a single open channel closes. This assumes that the average number of channels opening/closing in Δt is small. Then the probability that a single channel both opens and recloses (or vice versa) within a single time step can be ignored. To ensure this is the case, a suitably small Δt is required. More exact simulation methods are then not required [32]. Assuming that channels open and close independently of each other, the change in the number of open channels (and hence the value of f) can then be simulated by drawing from a binomial distribution. See S1 Text for details.

The number of each channel type changes from cell to cell. We assume that, for each channel type, there is some approximately constant area density so that the channel number scales with the square of the cell size. (In the following we always take cell size to refer to the diameter.) Based on measured single and total conductances, we estimate that in a typical cell of diameter 10 μm with g_BK = 0.5 nS there are about 5 BK channels, 200 SK channels, 640 K channels and 200 Ca channels (see Methods for details). We have checked that our conclusions are not altered if these values are different by up to a factor of two (see S1 Text).

Channel noise smooths the spiking-bursting transition

We are interested in excitable behaviour where the membrane potential increases significantly from its resting value, enters a depolarised state and later returns to the resting state. We call such behaviour an event and define its duration using a threshold for V (Fig 1A). Events are caused by interplay between the calcium and potassium currents. At the start of an event, Ca²⁺ ions flow into the cell and so increase V. Later, the slower potassium currents respond by flowing out of the cell, decreasing V and so terminating the event. Two broad categories of event are observed: spikes, where the potential returns quickly to its resting value without any oscillations, and bursts, where significant time is spent in the depolarised state, often with oscillations around some depolarised value. For concreteness we here define spikes as events with duration less than 100 ms and no oscillations. Events of either longer duration or with oscillations are classified as bursts. Other definitions are possible but are unlikely to affect our conclusions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Spikes, bursts and the effect of noise.

(A) The duration of an event is determined by the time that V spends above some threshold (chosen throughout as −45 mV). Events are classified as spikes if the event duration is less than 100 ms and there are no oscillations in V (e.g. first, yellow event). Events of longer duration or events with oscillations are identified as bursts (e.g. second, orange event). (B) The deterministic model with g_BK = 0.5 nS leads to pure periodic spiking. (C) The stochastic model with the same parameters shows that noise can convert some spikes into bursts. (D) Increasing g_BK to 1 nS leads to pure bursting in the deterministic model. (E) The equivalent stochastic model shows a mixture of spikes and bursts, with noise responsible for converting some bursts into spikes. (F) The deterministic model with g_Ca doubled to 4 nS simply sits at a depolarised steady state (g_BK = 0.5 nS). (G) However, the stochastic model with the same parameters shows bursts as noise continually drives the system away from the steady state.

https://doi.org/10.1371/journal.pcbi.1007769.g001

The deterministic model has both periodic spiking and bursting solutions [16]. For example, consider varying the total BK conductance g_BK which, for fixed single channel conductance g_1,BK, is equivalent to varying the number of BK channels N_BK. At low g_BK (only a few BK channels) we find purely spiking behaviour just as in the classic Hodgkin-Huxley model (Fig 1B). However, at large g_BK (many BK channels) all spikes disappear giving rise to pure bursting (Fig 1D). This is caused by the relatively fast response of BK channels (τ_BK < τ_n), which limits the initial membrane depolarisation. As a consequence, K channel activation is delayed, leading to longer events and oscillations in V around some depolarised value.

Turning on noise can dramatically change this behaviour. Although a simple implementation of noise has been studied before [16], this is the first time, as far as we are aware, that realistic channel noise with realistic numbers of channels has been considered in endocrine cells. For example, some of the spikes seen when g_BK is small can be converted into bursts (Fig 1C). Similarly, at large g_BK, spikes can be seen in addition to bursts (Fig 1E). Further, even cases where deterministically there are no events at all (Fig 1F) can acquire events when noise is turned on (Fig 1G).

Solutions with both spikes and bursts lead to the concept of the bursting fraction (BF), defined as the fraction of events over a long time period that are bursts. Much of this paper will be concerned with how BF depends on model parameters and cell size. Almost all deterministic solutions that have events correspond to either pure spiking (BF = 0) or pure bursting (BF = 1) so that, as the BK conductance g_BK gradually increases, there is a critical point where BF suddenly switches from 0 to 1 (Fig 2). The inclusion of noise smooths out this behaviour and leads to two important differences. Firstly, bursts start at a lower value of g_BK due to noise converting some spikes into bursts. Secondly, BF = 1 is never reached, even at very high g_BK: noise always converts some bursts to spikes. (This is in contrast to previous work where noise was treated less realistically, resulting in BF tending to 1 at large g_BK [16]). At first glance this seems counterintuitive: large g_BK corresponds to large numbers of BK channels and so should tend to the deterministic behaviour (BF = 1). However, although increasing g_BK increases N_BK, the single channel conductance g_1,BK is unchanged. Thus, the change in current due to one BK channel opening or closing is always g_1,BK(V − V_K), independent of N_BK, so that the proportion of bursts converted into spikes is similar for all large values of N_BK.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The bursting fraction against total BK conductance.

In the deterministic model (dashed curve), the bursting fraction (BF) shows a sharp jump from pure spiking (BF = 0) to pure bursting (BF = 1) at a critical value of g_BK. However, with channel noise (solid curve), this step function is smoothed out and no longer reaches pure bursting even for large g_BK. Here, g_1,BK has been fixed throughout so that g_BK is proportional to N_BK. Note that N_BK values only have a physical meaning in the stochastic model.

https://doi.org/10.1371/journal.pcbi.1007769.g002

The number of channels and the BK time constant

Here we first consider the effect of the total number of channels on the bursting fraction. To do this, we scale the numbers of all four channel types by a factor σ so that {N_Ca, N_K, N_SK, N_BK} → {σN_Ca, σN_K, σN_SK, σN_BK} whilst keeping the total conductances fixed and ensuring that σN_X is an integer. This is equivalent to scaling all four single channel conductances by 1/σ. Scaling of single channel conductances whilst keeping the total conductances fixed is unlikely to be biological but shows the effect of the noise level on electrical activity. As σ → ∞ we expect to recover the deterministic model.

First, consider the case when g_BK = 0.5 nS, which corresponds to pure spiking in the deterministic model. For the smallest values of σ, when there are only one or two BK channels, almost exactly half of events are spikes and half are bursts (Fig 3A). As σ increases, the bursting fraction gradually decreases to its deterministic value of BF = 0. Next, consider larger values of g_BK, which lead to pure bursting in the deterministic case. Again, for small σ, about half of events are spikes/bursts (Fig 3A). As σ increases, BF increases towards 1 as expected, but not before first initially dipping slightly. This means that, for deterministically-bursting values of g_BK, the least bursting corresponds to an intermediate number of channels, which in turn suggests that one way of controlling the bursting fraction involves careful choice of both the channel numbers and the single channel conductances.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The effect on bursting of the total number of channels and the BK time constant.

(A) For g_BK = 0.5 nS, corresponding to pure spiking in the deterministic model, BF gradually decreases as the total number of channels increases. However, for g_BK = 0.6 nS, corresponding to pure bursting deterministically, increasing σ leads to an initially drop in BF before an increase towards BF = 1. (B) Higher values of τ_BK lead to less bursting as BK channels start to act more like K channels. Conversely smaller τ_BK, i.e. even faster acting BK channels, lead to more bursting. (C) The overall sigmoidal BF-against-g_BK behaviour seen in Fig 2 is the same for other values of τ_BK. However, sufficiently large values of τ_BK only rarely lead to bursting even for very large total BK conductances.

https://doi.org/10.1371/journal.pcbi.1007769.g003

In the deterministic model, and experimentally, we have shown that the BK channel needs to have fast opening and closing kinetics in order to have a burst-promoting effect [9, 16]. We now examine the effect of the BK time constant, τ_BK in our stochastic model. Larger τ_BK values mean that the BK channels start to behave more like the slower-acting K channels. Thus they are less likely to convert spikes to bursts, resulting in lower BF (Fig 3B). Conversely, smaller τ_BK correspond to very fast-acting BK channels and a tendency for almost all events to be bursts. This is in agreement with the deterministic model [16]. The effect on the sigmoidal BF-against-g_BK behaviour in Fig 2 is similar: small τ_BK result in effectively pure bursting even for modest values of g_BK, whereas large τ_BK require larger g_BK before any noticeable bursting and, even for very large g_BK, only occasionally generate bursts (Fig 3C).

Individual channel types: BK versus other channels

The above results cover the realistic case of stochastic opening and closing in all four channel types that we consider here. However, with our model, we are able to determine the relative effect of noise in different channels. For example, we can examine whether noise in BK channels is more likely to cause bursting than noise in other channels. We do this by allowing some channel types to be described stochastically (as described above) and others deterministically (by considering the finite difference version of Eq (3)). In particular, we focus on comparing noise only in BK channels (with deterministic Ca, K and SK channels) with noise only in Ca, K and SK channels (with BK deterministic). We refer to the first case as “BK noise” and the second as “non-BK noise”.

First, consider the plot of the time derivative of the membrane potential () against the potential itself (V). For the case of pure spiking in the deterministic model (e.g. with g_BK = 0.5 nS) the -V trajectory is a single, non-intersecting loop (the dashed line in Fig 4A). This contrasts with the pure bursting case (e.g. with g_BK = 0.6 nS) where the trajectory self-intersects as V oscillates around a depolarised state. With the inclusion of stochastic channels, these trajectories become noisy. However, the precise nature of how noise affects the trajectories depends on exactly which channels are stochastic. For non-BK noise, the trajectory shows considerable fluctuation at all points during the cycle (Fig 4B). However, for BK noise, the trajectory is much smoother, with most of the fluctuation occurring while in the depolarised state (Fig 4A). This is related both to the fact that (V − V_K) is small when in the resting state and to the fact that, during most of the cycle, when V is close to the resting potential, there are rarely any BK channels open; this, in turn, is because the steady state activation of BK channels (f_∞ in Eq (4)) is very steep compared to that for other channels (since s_f is small compared to s_m and s_n). Only once an event has started and V begins to rise, does the stochastic nature of BK channels become apparent and lead to noticeable differences from the deterministic model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Relative noise contributions for BK versus non-BK channels.

(A) Plot of dV/dt against V with g_BK = 0.5 nS and noise only in BK channels (Ca, K and SK behave deterministically) for a short 600 ms period during which two events occur (a spike followed by a burst). The dashed line shows the entirely deterministic solution (with no noise in any channels), which corresponds to continuous spiking. The black circle labels the initial starting configuration. (B) As in (A) but with noise only in the Ca, K and SK channels (BK acts deterministically). (C) The BF-against-g_BK behaviour for different sources of noise showing that, whereas non-BK noise slightly smooths out the deterministic step function, the majority of the difference between the deterministic and fully stochastic models arises from noise in the BK channels.

https://doi.org/10.1371/journal.pcbi.1007769.g004

Now consider the maximum potential V_max during an event, which, due to noise, varies from event to event. For non-BK noise, V_max is similar for spikes and bursts. For example, for g_BK = 0.5 nS, the average of the V_max distribution is −5.6 mV for spikes and −5.9 mV for bursts, with a standard deviation of 1.1 mV in both cases. This should be compared to the deterministic result (corresponding to pure spiking), where V_max = −5.9 mV. The situation is different for BK noise, where the average V_max for bursts (−7.3 mV) is significantly lower than that for spikes (−4.8 mV). Again, the standard deviation (2.4 mV) is similar in each case. This difference between BK and non-BK noise also manifests itself in the BF-against-V_max histogram, which is approximately uniform for noise in non-BK channels and monotonically decreases towards 0 for BK noise. These results are similar for other values of g_BK and show that the way noise converts an event that would deterministically be a spike into a burst is different for BK and non-BK channels. For the case of BK noise, a burst is often caused by the sudden noisy opening of a BK channel as V is rising towards its maximum. This reduces V_max, which in turn delays K channel activation and so prolongs the event. Conversely, for non-BK noise, bursts are typically due to the noisy opening or closing of a non-BK channel only once V has reached its maximum and has started to drop back towards its resting value.

The difference between BK and non-BK noise can also be seen in the distribution of event durations. In all cases this distribution ranges from approximately 50 ms for the shortest spikes to about 250 ms for the longest bursts. For concreteness, consider g_BK = 0.5 nS, which shows both spiking and bursting events. For the case of BK noise, the distribution of event durations is unimodal and fairly uniform. However, for non-BK noise this distribution is distinctly bimodal, with a clear separation between events that are spikes and events that are bursts.

Finally, we examine the BF-versus-N_BK curve. The effect of either BK or non-BK noise lies somewhere between the deterministic result (dashed line in Fig 2) and the full noise case (solid line in Fig 2), with both smoothing the deterministic step function to a greater or lesser extent (Fig 4C). However, non-BK noise only slightly differs from the deterministic case (around the spiking-bursting transition), with BF tending towards pure bursting for large N_BK. In contrast, the case of BK noise looks almost identical to the full noise case. Thus, although noise in non-BK channels accounts for most of the fluctuation in the -V trajectory (Fig 4B), it is BK channel stochasticity that determines the fraction of events that are bursts and typically causes events that would be spikes (or equivalently bursts) in the deterministic model to sometimes instead become bursts (or spikes). As expected, the fact that BK channels have a greater effect than other channels is due to their scarcity and relatively large single channel conductance.

The role of cell size

Mathematical models of endocrine pituitary cells often neglect the role of cell size, even though there are substantial variations in the size of these cells, with cell diameters ranging from 7 to 15 μm [33, 34]. Cell size varies not only between cell types but also within a given cell type and, for many organisms, depends on age, sex and phase within the reproductive cycle [35–37]. Further, pituitary adenomas and other pituitary tumours are often associated with changes in cell size [38].

We now explore the effect of cell size on electrical activity. We assume that cell size can be described by some effective radius R. It is not necessary that the cell is spherical, only that its membrane area and volume are proportional to R² and R³ respectively. As our base case, we take a cell with diameter 2R = 10 μm and consider scaling the radius by a factor λ so that R → λR. Then the membrane area and cell volume scale by λ² and λ³ respectively.

This scaling is likely to affect only some parameters. For example, we assume that the reversal potential and current gating parameters (voltage and calcium midpoints and slopes) are unaffected. However, other parameters, such as the conductances and membrane capacitance, will change. In particular, assuming a constant area density of membrane channels, we expect that both the numbers of channels N_X and therefore the conductances g_X scale with the membrane area. This also applies to the leak conductance. The membrane capacitance, C, which is the charge that the membrane stores per unit electric potential, also scales with membrane area.

Since the charge-to-molar conversion in the [Ca] equation, α, converts a current to an intracellular Ca²⁺ concentration, it follows that α scales inversely with the cell volume. In addition, assuming calcium extrusion occurs through the Na⁺/Ca²⁺ exchanger (NCX) and the plasma membrane calcium ATPase (PMCA), both of which operate through the membrane, we expect the Ca²⁺ extrusion rate k_c to scale like the membrane area over the cell volume. This is because the extruded calcium leaves through the membrane but affects the calcium concentration throughout the whole cell. (We have here assumed that the calcium dynamics occur throughout the whole cell volume. Although we believe this to be the most likely scenario, it is possible instead that the relevant changes in [Ca] are confined to a thin shell region adjacent to the membrane. This would lead to a different, although still nontrivial, dependence on cell size. See S1 Text for a full discussion of this.)

Thus, in summary, scaling R → λR implies that (9) Since the membrane capacitance and conductances both scale with λ², the equation for the membrane potential (Eq (1)) is unaffected by the R → λR scaling. This is also the case for the m, n, s and f gating variables. However, the [Ca] equation (Eq (5)) acquires an overall 1/λ scaling on the right-hand side. It is worth noting that, if we only focus on the six model variables (V, m, n, s, f, [Ca]) and ignore other quantities such as the currents, this can equivalently be seen as a 1/λ scaling to the f_c parameter (the fraction of intracellular calcium present as free ions). Then the overall effect of scaling R → λR is equivalent to simply scaling (10) Of course, this does not mean that the intracellular ratio of free to bound calcium depends on the cell size, simply that this scaling of f_c gives an equivalent mathematical model.

Cell size in the deterministic model

We now study the consequences of cell size on spiking and bursting behaviour. First, we investigate the deterministic model. Since this is equivalent to the limit N_X → ∞, the only effect of scaling R → λR is to send f_c → f_c/λ. This sharply contrasts with the classic deterministic Hodgkin-Huxley model, where the absence of an equation for [Ca] means there is no effect of scaling R and so identical dynamics for all cell sizes. It is precisely because this is not the case in the present model that makes the role of cell size, even in the deterministic version, an interesting topic.

For concreteness, consider the case where g_BK = 0.5 nS at λ = 1, which produces spiking in the deterministic model. First, we examine what happens as the cell size is increased (λ > 1), with a corresponding increase in g_BK (Eq (9)). This is equivalent to reducing f_c and so effectively results in slower calcium dynamics. As a consequence, once an event starts, the intracellular calcium concentration increases slower than for λ = 1, which causes s_∞ (and so s) to also increase slower. In turn, this means that I_SK rises more slowly, so that V remains longer in the depolarised state and so leads to longer events. Modest increases in λ still correspond to spiking (although with longer event durations), but for sufficiently large λ there is a sudden transition to pure bursting. Further increasing λ results in longer and longer bursts. It is noteworthy that the mechanism by which spiking transitions to bursting when λ is increased is different to that when g_BK is increased by itself: the former is related to slower calcium dynamics, whereas the latter involves a larger BK current and so lower V_max.

Now consider the case of smaller cells (λ < 1). This is equivalent to quicker calcium dynamics, which in turn leads to shorter events. This means that spiking still continues for smaller cells, although with shorter durations. Only for λ < 0.02 (which corresponds to unnaturally small, submicron-sized cells) does the spiking behaviour break down, with the cell instead entering a permanently depolarised state with V ≈ −45 mV. These results for both smaller and larger cells are summarised in Fig 5A, which shows that the transition to bursting is accompanied by a sharp increase in event duration. The non-monotonic behaviour within the bursting region is due to oscillations of V within the depolarised state: each maximum corresponds to the appearance of an extra oscillation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. The effect of cell size in the deterministic model.

(A) Event duration against cell size for the case where g_BK = 0.5 nS at λ = 1. The spiking behaviour switches sharply to bursting at around λ = 1.35. Light shading: pure spiking; dark shading: pure bursting. (B) The critical BK conductance for switching between spiking and bursting as a function of the cell size. The blue curve (left axis) shows the total BK conductance within the cell (g_BK), whereas the red curve (right axis) shows the BK conductance area density (g_BK/4πR²).

https://doi.org/10.1371/journal.pcbi.1007769.g005

Finally, we consider how, in the deterministic model, the transition between spiking and bursting depends on cell size. For a given λ, there is a critical conductance, , at which spiking behaviour sharply switches to bursting (see the dashed curve in Fig 2). Since the total BK conductance scales with cell size (Eq (9)), it makes more sense to consider the g_BK area density (i.e. the BK conductance due to a unit area of membrane), defined as g_BK/4πR². To calculate this, we identify λ = 1 with a diameter of 10 μm, although our qualitative results are not dependent on this choice. In Fig 5B we plot how both the critical g_BK and the critical g_BK density depend on λ. As expected, larger cells (with more channels) require larger to transition from spiking to bursting. However, interestingly, the critical g_BK density (i.e. ) decreases with increasing cell size, suggesting that larger cells are more likely to exhibit bursting than smaller cells. This can again be understood from the effectively slower calcium dynamics in larger cells, which leads to longer events and means that the spiking-to-bursting transition occurs at lower g_BK densities.

Cell size in the stochastic model

We now consider the effect of cell size in the more realistic stochastic model. Only values of λ that lead (without rounding) to an integer number of channels (for all channel types) are considered. In particular, since we take N_BK = 5 at λ = 1, we do not consider cells smaller than λ² = 1/5; such cells would not even have one BK channel. The effect of changing λ now not only affects the rate of calcium dynamics, but also changes the number of channels (Eq (10)). In Fig 6A we plot example voltage traces for five different cell sizes. As usual, we set g_BK = 0.5 nS at λ = 1. As expected, there is less noise in the V time course for larger cells due to the greater number of channels. Further, larger cells produce more bursts and generally longer events.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. The effect of cell size in the stochastic model.

(A) Examples of the membrane potential from simulations of five different cell sizes, showing mostly spiking for the smallest cells (λ = 0.45) and gradually becoming almost complete bursting for large cells (λ = 3.16). For each trace, the y-axis ranges from −70 to 0 mV. (B) The bursting fraction as a function of the cell size based on a 10 μm cell (λ = 1) having 5 BK channels. The dashed line shows the equivalent plot for the deterministic model. (C) The bursting fraction as a function of the g_BK area density for three cell sizes. Panel (B) corresponds to BF measured along the dashed vertical line in Panel (C).

https://doi.org/10.1371/journal.pcbi.1007769.g006

To examine this further, we plot the bursting fraction against cell size (Fig 6B). It is worth noting that, based on a 10 μm cell for λ = 1, the extremes of this plot correspond to unrealistically small (<5 μm) and large (>30 μm) cells. As cell size increases, BF increases towards pure bursting. This is for exactly the same reasons as in the deterministic model: increasing λ corresponds to slower calcium dynamics, which in turn leads to a slower rise in the SK current, giving rise to longer events that are more likely to be bursts. Similarly, up to a point, smaller cells result in shorter events and so lower bursting fractions. However, interestingly, for the smallest cells, BF again starts to rise. This is because smaller cells have fewer channels, meaning that stochastic effects become more important. BF is then determined by competing effects between the calcium dynamics and the number of channels: for small cells, in addition to slow calcium dynamics, which tends to reduce event duration (and so reduce BF), there are also fewer channels (in particular BK channels), so that channel noise is more likely to convert would-be spikes into bursts (and so increase BF). This increase in bursting for the smallest cells is completely missed in the deterministic model and shows that the most spiking (smallest BF) is found in intermediate-sized cells.

We also plot BF against the g_BK density for various cell sizes (Fig 6C). Fig 6B can be thought of as Fig 6C along the dotted vertical line. In each case, we obtain the familiar smoothed Heaviside function with, except for the smallest g_BK densities, more bursting in larger cells. It is notable that the g_BK density at which bursts first appear (i.e. the point where BF first differs noticeably from 0) is approximately the same for all three cell sizes. In addition to our model predicting that large cells burst more, we also find that, in the smallest cells, at most approximately half the events can be bursts even for the highest g_BK densities, in stark contrast to the deterministic model.

The effect of opening and closing individual channels

Ion channel noise is able to convert would-be spikes into bursts and vice versa. In this section we examine, both in the deterministic and stochastic models, how this is related to the opening and closing of individual ion channels. We determine at what point during an action potential noise is most effective in converting a would-be spike into a burst and vice versa. For concreteness we focus on the case g_BK = 0.5 nS, which corresponds to pure spiking in the deterministic model.

In particular, we study the effect of manually perturbing the system at a particular point during an event and in a particular channel. The perturbation involves either opening/closing a given number of channels or injecting a current for a short fixed period of time. The point to apply the perturbation is determined from a potential threshold V*, with the perturbation applied the first time during each event that the potential crosses V*. In turn V* can be mapped, using the deterministic spike action potential, onto a time before or after V_max. For negative times (i.e. times whilst the deterministic potential is rising towards V_max), the perturbation is applied the first time V rises above V*. Conversely, for positive times (i.e. times whilst the deterministic potential is dropping from V_max), the perturbation is applied the first time (since V_max) that the potential drops below V*. This ensures that there is at most one perturbation per event.

In the stochastic model, for a given channel type, the size of the fluctuations of the number of open channels depends on the total number of channels and so on the channel type. For example, consider the BK channels. At any given time, at least if we assume a quasi-steady state, we expect the number of open BK channels to be approximately binomially distributed with mean f_∞N_BK and variance f_∞(1 − f_∞)N_BK (see S1 Text). Thus the maximum standard deviation (which occurs when f_∞ = 0.5) is . This suggests that, to make the size of the perturbation comparable across different channel types, we should take the perturbation for channel X to be . In particular, this means that we take a perturbation to consist of simultaneously opening or closing 7 Ca channels, 13 K channels, 7 SK channels or 1 BK channel.

Consider first the deterministic model. Based on Eq (2), we implement a perturbation in channel X by injecting a current g_X(q/N_X)(V − V_X) for 5 ms. Here q is the change in the number of open channels: positive/negative q corresponds to opening/closing q channels. Whether the perturbation transforms spikes into bursts depends both on q and on the time the perturbation is applied (i.e. on V*). In Fig 7A we plot the regions where bursting occurs along with the deterministic spike potential.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. The effect of opening and closing individual channels.

(A) Ranges over which opening or closing channels in the deterministic model causes spikes to become bursts. The blue curve shows the (non-perturbed) deterministic V profile, which corresponds to a spike. The horizontal axis shows the point during the action potential when the perturbation (here implemented as an extra current) is applied. Red/yellow regions: times from V_max when adding/removing a current converts spiking to bursting. Top/bottom: perturbation applied in BK/Ca channels. In each case the perturbation is applied for 5 ms. (B) Effect on BF in the stochastic model of opening/closing a single BK channel at different points during the action potential. The horizontal axis shows the point when the perturbation is applied: a time from V_max corresponds to a particular value of the potential in the deterministic model (blue curves in (A)); the perturbation is applied when the potential first rises above or drops below this value. Blue horizontal dashed line: value of BF with no perturbation.

https://doi.org/10.1371/journal.pcbi.1007769.g007

For potassium channels (upper panel in Fig 7A), there is a region around V_max where transiently opening channels leads to bursting. This works in much the same way that BK channels cause bursting: the perturbation lowers V, which delays subsequent activation of the repolarising currents, and so extends the event duration. This effect only applies within 5 to 10 ms of V_max: perturbations that are applied too early or too late make little difference to the event profile. Similarly there is a region after V_max (as V is dropping back towards the resting state) when closing potassium channels leads to bursts. In this case the applied current acts to increase V and so maintains the cell in the depolarised state for longer. The size of the regions that lead to bursting varies depending on the channel type. Opening or closing a single BK channel leads to wider regions than opening or closing seven SK channels, which in turn leads to wider regions than opening or closing 13 K channels. This difference boils down to the relative size of g_X(q/N_X), which is proportional to .

For calcium channels (lower panel in Fig 7A), we instead find the opposite effect: manually opening Ca channels leads to bursting only for a time after V_max, whereas closing Ca channels causes bursts when V is around its maximum. This corresponds to the fact that, for V near V_max, I_Ca is negative (calcium ions flow into the cell) whereas I_K/SK/BK is positive (potassium ions flow out of the cell). Interestingly the range within which opening Ca channels after V_max causes bursting is significantly larger than that for closing BK channels.

We now consider perturbations in the full stochastic model and determine how these affect BF as a function of when they are applied. In particular, we plot BF against the time in the deterministic model when the potential first crosses V* (see Fig 7B for the BK case). Rather than using currents, perturbations are now implemented as the direct opening or closing of a given number of channels. When potassium channels are opened/closed before or around V_max, BF increases/decreases, whereas when they are opened/closed after V_max, BF is decreased/increased. Again, the opposite is seen for calcium channels. The point at which perturbations make no change to BF depends on channel type and does not correspond to V_max. For example, for BK channels, a perturbation applied about 8 ms after V_max does not affect BF. This does not mean that these perturbations do not cause spikes to become bursts and vice versa, simply that the number of spikes converted to bursts is offset by the number of bursts converted to spikes.

Not surprisingly, the size of perturbation required for an observable change to BF depends on the channel type: perturbations of only one K or BK channel notably alter BF, whereas more than 40 Ca or SK channels must be simultaneously opened or closed for any noticeable effect. This is related to the very short time constants for Ca and SK channels (τ_m and τ_s), which mean that modest perturbations in these channels are quickly rectified. This is not the case for the K and BK channels, where τ_n and τ_BK are much longer.

The size of the shift in BF often depends on whether the perturbation is before or after V_max. For example, for K and Ca channels, perturbations have the largest effect when applied after V_max, with a much smaller effect for perturbations before V_max. This is not the case for BK channels with approximately equal shift in BF (see Fig 7B). The situation for SK channels is difficult to determine since even simultaneously opening or closing 100 channels has little effect on BF. These results match well with those in the “Individual channel types: BK versus other channels” section where we considered the effect due to noise in only certain channel types at a time.

Finally, it is worth noting the small local maximum when the perturbation involves closing a BK channel at V_max (yellow curve in Fig 7B). The equivalent is not seen when a BK channel is opened. This is because, for V* near V_max, not all events reach a high enough V for the perturbation to be applied. In particular, when a BK channel is closed near V_max the effect is to convert some bursts into spikes and so reduce BF. However, because bursts tend to have lower V_max than spikes, some bursts avoid the perturbation and so remain bursts, meaning that BF is not reduced as much as might initially be thought. The opposite effect due to opening a BK channel near V_max is not seen since would-be spikes are more likely to cross the relatively high threshold value that triggers the perturbation.

In summary, we have been able to shed some light on the role channel noise plays in determining whether an event is a spike or a burst, and how this depends on channel type and fluctuation size. As expected, BK channels have the most significant effect, with the opening or closing of even a single channel often substantially affecting BF. In particular, the greatest effect occurs for BK perturbations either about 4 ms before V_max or 15 ms after. By contrast ten or more K channels must simultaneously be opened to observe a similar effect, and even then only for perturbations applied after V_max. Finally, the short time constants for the Ca and SK channels mean that noise in these channels typically has little effect on whether an event is a spike or a burst.

Robustness of results to parameter variations

The model parameters are likely to vary from cell to cell. This will be the case even for a particular cell type and when cell size is taken into account. To determine whether this affects our results, we used a GPU to run our simulations for a range of parameter sets [39]. This allowed us to efficiently check hundreds of thousands of different parameter choices. In particular, we chose to vary the four single channel conductances (g_1,Ca, g_1,K, g_1,SK, g_1,BK), the leak conductance and reversal potential (g_l, V_l), and the calcium extrusion rate k_c. For each simulation, we picked each parameter randomly from a uniform distribution with a range ±50% of the parameter’s default value.

Each parameter set can then be classified into one of four basic types—depolarised, hyperpolarised, noisy steady state or event-containing—based on the behaviour of the membrane potential V (once the initial transient period has passed). The depolarised and hyperpolarised cases are defined by the size of variations in V being less than 10 mV, with depolarised (hyperpolarised) corresponding to whether the mid-range of V is above (below) −50 mV. Noisy steady states encompass a number of different behaviours: (i) cases with no events, (ii) cases where the range of V lies between 10 mV and 35 mV, and (iii) cases where the average event duration is much longer than the average inter-event time. All other cases are labelled as event-containing. We further classify the event-containing cases as either pure spiking (BF = 0), almost pure spiking (0<BF<0.05), pure bursting (BF = 1), almost pure bursting (0.95<BF<1) or mixed (0.05≤BF≤0.95). See S1 Text for full details.

Based on 700,000 parameter sets, we find that 92.4% contain events, with only 0.1% depolarised, 3.0% hyperpolarised and 4.5% noisy steady states. Of those that contain events, 3.9% correspond to pure spiking, 5.8% to almost pure spiking, 6.0% to pure bursting, 6.8% to almost pure bursting, and 77.5% to mixed cases.

We now examine the sigmoidal behaviour of BF against N_BK. Above, we found that BF increases with N_BK (with g_1,BK fixed) but typically saturates before reaching pure bursting (Fig 2). We can now test whether this is also the case for other parameter choices. For each randomly-chosen parameter set, we consider a range of values for N_BK (with a corresponding change in g_BK). We only study cases where all choices of N_BK are event-containing and where there is a large enough BF range to make comparisons meaningful. For example, cases where nearly all N_BK values corresponded to pure spiking were not included. Of the selected parameter sets, over 98% (141,449 of 143,870) exhibit the sigmoidal behaviour shown in Fig 2. Note that we allow small decreases in BF for large N_BK as happens in Fig 2. Further, we find that the large N_BK limit of BF is usually less than unity, with an average of 0.80 and a standard deviation of 0.12. Thus, for most parameter choices, even very large numbers of BK channels do not lead to pure bursting.

The few cases where BF does not increase with N_BK are interesting to study, and there are even cases where BF decreases with increasing N_BK. For example, this occurs if the maximal conductance for Ca channels is increased to g_Ca = 3.5 nS and that for K channels decreased to g_K = 1.8 nS. In the deterministic model, this parameter choice produces a depolarised steady state for N_BK < 18 and pure bursting for larger N_BK. However, the inclusion of channel noise causes small N_BK to now correspond to pure bursting, with BF gradually dropping to around 0.4 for large N_BK.

Finally, we examine our result that larger cells exhibit more bursting. Again we choose random parameter sets, but now compare three cell sizes (λ = 0.5, 1.0 and 1.5) at two different channel densities (corresponding to N_BK = 4 and N_BK = 20 when λ = 1.0). As before, we neglect cases where not all simulations contain events. Of the remaining cases, 89.8% (144,670 of 161,116 parameter sets) for the first channel density and 99.7% (153,781 of 154,272) for the second show BF increasing with cell size.

Parameter values

Unless otherwise stated, we use the following parameter values [16]. C = 10 pF, g_Ca = 2 nS, g_K = 3.2 nS, g_SK = 2 nS, g_BK = 0.5 nS, g_l = 0.2 nS, V_Ca = 60 mV, V_K = −75 mV, V_l = −50 mV, τ_m = 0.1 ms, τ_n = 30 ms, τ_s = 0.1 ms, τ_BK = 5 ms, v_m = −20 mV, s_m = 12 mV, v_n = −5 mV, s_n = 10 mV, v_f = −20 mV, s_f = 2 mV, k_s = 0.4 μM, f_c = 0.01, α = 0.0015 μM/fC, k_c = 0.12 ms⁻¹.

Estimation of channel numbers

The number of channels is likely to depend on the cell size. For simplicity we assume that, for a given channel type, there is a constant number of channels per unit membrane area. In particular, for concreteness, we take a typical cell to have diameter 10 μm and contain 200 Ca channels, 640 K channels, 200 SK channels and 5 BK channels. These are based on literature values for the single channel conductances of g_1,Ca = 10 pS, g_1,K = 5 pS, g_1,SK = 10 pS and g_1,BK = 100 pS [49, 50]. We have checked that our conclusions are unchanged if these values are decreased by up to 50% or increased by up to 100% (see S1 Text).

Numerical simulations

The system was solved numerically using the Euler method with step size Δt = 0.01 ms. Each time step involved (i) calculating the currents (Eq (2)), (ii) applying any perturbation, (iii) updating V and [Ca] (Eqs (1) and (5)), (iv) calculating the steady-state activation functions (Eq (4)), and (v) stochastically determining how many channels open and close during the time step. We checked (for a range of parameter values) that smaller values of Δt, down to Δt = 1 × 10⁻⁴ ms, do not lead to noticeably different results. For full details, see S1 Text.