Mathematical Modeling of Heterogeneous Electrophysiological Responses in Human β-Cells

Electrical activity plays a pivotal role in glucose-stimulated insulin secretion from pancreatic -cells. Recent findings have shown that the electrophysiological characteristics of human -cells differ from their rodent counterparts. We show that the electrophysiological responses in human -cells to a range of ion channels antagonists are heterogeneous. In some cells, inhibition of small-conductance potassium currents has no effect on action potential firing, while it increases the firing frequency dramatically in other cells. Sodium channel block can sometimes reduce action potential amplitude, sometimes abolish electrical activity, and in some cells even change spiking electrical activity to rapid bursting. We show that, in contrast to L-type -channels, P/Q-type -currents are not necessary for action potential generation, and, surprisingly, a P/Q-type -channel antagonist even accelerates action potential firing. By including SK-channels and dynamics in a previous mathematical model of electrical activity in human -cells, we investigate the heterogeneous and nonintuitive electrophysiological responses to ion channel antagonists, and use our findings to obtain insight in previously published insulin secretion measurements. Using our model we also study paracrine signals, and simulate slow oscillations by adding a glycolytic oscillatory component to the electrophysiological model. The heterogenous electrophysiological responses in human -cells must be taken into account for a deeper understanding of the mechanisms underlying insulin secretion in health and disease, and as shown here, the interdisciplinary combination of experiments and modeling increases our understanding of human -cell physiology.


Introduction
Glucose-stimulated insulin secretion from human pancreatic bcells relies on the same major signaling cascade as their rodent counterparts, with electrical activity playing a pivotal role. Following metabolism of the sugar, ATP-sensitive potassium channels (K(ATP)-channels) close in response to the elevated ATP/ADP-ratio, which triggers action potential firing and Ca 2z -influx through voltage-gated calcium channels. The resulting increase in intracellular calcium leads to insulin release by Ca 2z -dependent exocytosis [1][2][3][4]. However, the electrophysiological properties of human and rodent b-cells show important differences, e.g., with respect to their palette of expressed Ca 2z -channels and the role of Na z -channels, which contribute to electrical activity in human but not in rodent b-cells [1,3].
Mathematical modeling has played important roles in studying the dynamics of electrical activity in rodent b-cells [5,6], and could plausibly aid in understanding the electrophysiological responses in human b-cells, and how they might differ from rodent cells.
Recently, the first model of electrical activity in human b-cells [7] was constructed from careful biophysical characterizations of ion channels in human b-cells, mainly from Braun et al. [3]. The model [7] included Na z -channels, three types of Ca 2z -channels, an unspecified leak-current, and several K z -channels: delayed rectifier (Kv) K z -channels, large-conductance (BK) Ca 2z -sensitive K z -channels, human ether-a-go-go (HERG) K z -channels as well as K(ATP)-channels. Recently evidence for small conductance (SK) Ca 2z -sensitive K z -channels in human b-cells was published [4,8], a current not included in the mathematical model [7].
The model [7] was shown to reproduce, depending on parameter values, spiking or rapid bursting electrical activity, which could be modified in accordance with a series of experiments by simulating pharmacological interventions such as ion channel blocking. These experiments were in general straightforward to interpret, also without a model. For example, the facts that blocking depolarizing Na z -or Ca 2z -currents slowed or abolished electrical activity [3] are as one would expect.
Here, we extend the previous model for human b-cells [7] by including SK-channels and Ca 2z dynamics, and show that the model now has reached a level of maturity that allows us to get insight in less intuitive experimental findings. We find experimentally that SK-channels in some cells play an important role in controlling electrical activity, while they have virtually no effect in other cells. Using the extended version of the model, we show that this difference can be explained by differences in the excitability of the cells. Moreover we find that SK-channels can substitute for HERG-channels in controlling rapid bursting. We also show that blocking Na z -channels in some cells can transform spiking behavior into rapid bursting, in contrast to the usual effect of Na z -channel blockers, which in general reduce or abolish spiking behavior [3,9,10]. Using our model we suggest that this happens in cells with a large Na z -current and that BK-channels play a prominent role. In addition, we suggest that SK-channels might underlie the surprising result that blocking depolarizing P/Q-type Ca 2z -channels enhances electrical activity, in contrast to the effect of L-or T-type Ca 2z -channel antagonists, which reduce excitability and electrical activity [3]. Our model is then used to investigate paracrine effects of c-aminobutyric acid (GABA) and muscarinic signaling on electrical activity. Finally, we show experimentally slow oscillations in electrical activity that might underlie pulsatile insulin secretion from human pancreatic islets, and by adding an oscillatory glycolytic component [11] to the electrophysiological model, we simulate such slow bursting patterns.

Results
To investigate a series of experimental observations, we have extended our previous model of electric activity in human b-cells [7] by including several additional components of human b-cell physiology, as described in the following, and in greater details in the Methods section. The mathematical modeling was carefully based on experimental data, as was the development of the core electrophysiological part modeled previously [7]. The extended model includes small conductance Ca 2z -activated K z -channels (SK-channels), which are expressed in human b-cells [4,8]. The size of the SK-current was estimated from experimental measures [8]. We made a special effort to carefully model the submembrane dynamics of Ca 2z , since SK-channels are controlled by the submembrane Ca 2z concentration (½Ca 2z mem ), which reacts rapidly to each action potential so that activation of SK-channels might influence the generation and shape of action potentials during spiking electrical activity. In order to study paracrine signalling, our extended model also includes currents due to caminobutyric acid (GABA) activation of GABA A receptors, which are ligand-gated Cl { channels operating in human b-cells [12]. Finally, a glycolytic oscillator [11] has been added to the model to account for slow oscillations in ATP levels in human b-cells [13,14], which have been suggested to underlie slow patterns of electrical activity, Ca 2z oscillations, and pulsatile insulin release.
Summarizing, the new version of the model now includes components from glucose metabolism, additional electrophysiological components (SK-channels and GABA A receptors), and Ca 2z dynamics, leading to a global model of human b-cell physiology, which, importantly, is based as far as possible on published data from human b-cells.

SK channels
When stimulated by glucose, human b-cells show electrical activity [1,3]. Human b-cells express SK-channels [4,8], which might participate in controlling electrical activity. To study the role of SK-channels in human b-cells, we included SK-channels and Ca 2z dynamics in our previous model [7]. The new model with standard parameters produces spiking electrical activity (Fig. 1A), which is virtually unaffected by setting the SKconductance g SK~0 nS/pF simulating SK-channels block. This model prediction was confirmed by our experimental data, and was also observed in at least one cell by Jacobson et al. [8]. Fig. 1B shows an example of spiking electrical activity in a human b-cell stimulated by 6 mM glucose, where addition of the SK1-3 channel blocker UCL 1684 (0.2 mM) did not affect the spiking pattern. Unchanged or marginal effects on electrical activity were also seen with a specific SK4 channel antagonist, TRAM-34 (1 mM, Fig. 1C). However, in some cells TRAM-34 application increased the action potential dramatically (Fig. 1D) in agreement with observations with the SK-channel antagonist apamin [8]. Note that before SK-channel block, the cell in Fig. 1D was almost quiescent, and fired action potentials very infrequently and randomly. This increase in spike frequency can be simulated by a stochastic version of the model. By including noise in the K(ATP) current, an otherwise silent cell produces infrequent action potentials evoked by random perturbations (Fig. 1E). When the SK-conductance is set to 0 nS/pF, the cell starts rapid action potential firing driven by the underlying deterministic dynamics. The model analysis indicates that this mechanism only works if the cell is very near the threshold for electrical activity in the absence of the SK-channel antagonist. Düfer et al. [15] suggested a similar, important role for SK4 channels in promoting electrical activity in murine b-cells at subthreshold glucose concentrations. Summarizing, cell-to-cell heterogeneity can explain the differences seen in the electrophysiological responses to SK-channel antagonists.
In addition to spiking electrical activity, human b-cells often show rapid bursting, where clusters of a few action potentials (active phases) are separated by hyperpolarized silent phases [1,4,9,10,16] (Fig. 2A). The extended model presented here can also reproduce this behavior (Fig. 2B) as could the previous version of the model [7], where the alternations between silent and active phases were controlled by HERG-channels. In contrast, in the present version of the model the rapid burst pattern (Fig. 2B, upper trace) can be controlled by SK-channels, which in turn are regulated by ½Ca 2z mem and ultimately by bulk cytosolic Ca 2z

Author Summary
Insulin is a glucose-lowering hormone secreted from the pancreatic b-cells in response to raised plasma glucose levels, and it is now well-established that defective insulin secretion plays a pivotal role in the development of diabetes. The b-cells are electrically active, and use electrical activity to transduce an increase in glucose metabolism to calcium influx, which triggers insulin release. Experimental and theoretical studies on b-cells from rodents have provided valuable insight in their electrophysiology. However, human b-cells differ from their rodent counterparts in several aspects including their electrophysiological characteristics. We show that the electrophysiological responses in human b-cells to a range of experimental manipulations are heterogeneous. We extend a previous mathematical model of electrical activity in human b-cells to investigate such heterogeneous and nonintuitive electrophysiological responses, and use our findings to obtain insight in previously published insulin secretion measurements. By adding a glycolytic component to the electrophysiological model, we show that oscillations in glucose metabolism might underlie slow oscillations in electrical activity, calcium levels and insulin secretion observed experimentally. We conclude that the interdisciplinary combination of experiments and modeling increases our understanding of human b-cell physiology and provides new insight in b-cell heterogeneity.
levels (½Ca 2z c ). The simulated cytosolic Ca 2z concentration shows the characteristic sawtooth pattern (Fig. 2B, lower trace) of a slow variable underlying bursting [17,18]. Thus, as in the pioneering model by Chay and Keizer [19], ½Ca 2z c increases during the active phase and activates SK-channels, which eventually repolarize the cell. During the silent phase ½Ca 2z c decreases and SK-channels close, allowing another cycle to occur.

Na z channels
Blocking voltage-dependent Na z -channels in human b-cells showing spiking electrical activity with tetrodotoxin (TTX) typically reduces the action potential amplitude by ,10 mV, and broadens its duration [3,9,10] (Fig. 3A). The previous version of the model [7] could reproduce these results, though the reduction in peak voltage was slightly less than observed experimentally. The inclusion of SK-channels in the model leads to a greater reduction in the spike amplitude (Fig. 3B, upper trace) when Na z -channels are blocked. This improvement is because of a mechanism where the slower upstroke in the presence of Na zchannel blockers allows submembrane Ca 2z to build up earlier and to higher concentrations (Fig. 3B, lower trace), and consequently to activate more SK-channels, which in turn leads to an earlier repolarization reducing the action potential amplitude. In other experiments (Fig. 3C) [16], TTX application suppresses action potential firing. In agreement, simulated spiking electrical activity can be suppressed by TTX application if the cell is less excitable because of, for example, smaller Ca 2z -currents ( Fig. 3D, upper, black trace). Before TTX application, the simulated cell had less hyperpolarized inter-spike membrane potential (*{61 mV; Fig. 3D) compared to the simulation with default parameters (*{70 mV, Fig. 3B). This finding is in accordance with experimental recordings (compare Fig. 3A and 3C). The cessation of action potential firing leads to a reduction in simulated ½Ca 2z mem (Fig. 3D, lower, black trace). The model predicts that spiking, electrical activity can continue in presence of TTX even in less excitable cells, e.g., with lower depolarizing Ca 2z -currents, if the hyperpolarizing K(ATP)-current is sufficiently small (Fig. 3D, upper, gray trace). In this case, ½Ca 2z mem is nearly unchanged (Fig. 3D, lower, gray trace). Hence, it is the relative sizes of the depolarizing and hyperpolarizing currents that determine whether TTX application silences the cell or allows the cell to remain in a region where action potential firing continues. The model thus predicts that in some cells, which stop firing action potentials in the presence of TTX, increased glucose concentrations or sulfonylureas (K(ATP)-channel antagonists) could reintroduce spiking electrical activity.
More surprisingly, TTX application can change spiking electrical activity to rapid bursting in some cells (Fig. 3E). This behavior can also be captured by the model (Fig. 3F). To simulate this behavior it was necessary to increase the size of the Na zcurrent. Without TTX, the big Na z -current leads to large action potentials, which activate sufficient BK-current to send the membrane potential back to the hyperpolarized state, allowing a Figure 1. Heterogenous responses to SK-channel block. Note the differences in time-scales. A: Simulation, with default parameters, showing no effect of SK-channels block (g SK~0 nS/pF during the period indicated by the gray bar). B: Experimental recording of spiking electrical activity in the same human b-cell before (left) and during (right) application of the SK1-3 channel antagonist UCL-1684 (0.2 mM). C: Experimental recording of spiking electrical activity in the same human b-cell before (left) and during (right) application of the SK4 channel antagonist TRAM-34 (1 mM), which had little effect on the action potential frequency in this cell. D: Experimental recording of spiking electrical activity in the same human b-cell before (left) and during (right) application of the SK4 channel antagonist TRAM-34 (1 mM), which accelerated the action potential frequency in this cell. E: Stochastic simulation reproducing the dramatic effect of SK-channels block (g SK~0 nS/pF during the period indicated by the gray bar). Other parameters took default values, except g KATP~0 :0175 nS/pF. doi:10.1371/journal.pcbi.1003389.g001 new action potential to form. With Na z -channels blocked, there is insufficient depolarizing current to allow full action potentials to develop. In consequence, less BK-current is activated (Fig. 3F, lower trace), and the membrane potential enters a regime with more complex dynamics where smaller spikes appear in clusters from a plateau of ,240 mV. The change to bursting activity leads to a notable increase in simulated ½Ca 2z mem (Fig. 3F, middle trace).

Ca 2z channels
High-voltage activated L-and P/Q-type Ca 2z -currents are believed to be directly involved in exocytosis of secretory granules in human b-cells [1,3,4,20,21]. Blocking L-type Ca 2z -channels suppresses electrical activity [3], which is reproduced by the model (Fig. 4A) [7], and the lack of electrical activity is likely the main reason for the complete absence of glucose stimulated insulin secretion in the presence of L-type Ca 2z -channel blockers [3]. Thus, L-type Ca 2z -channels participate in the upstroke of action potentials and increases excitability of human b-cells.
In contrast, and surprisingly, application of the P/Q-type Ca 2z -channel antagonist v-agatoxin IVA does not block or slow down electrical activity, but leads to an increased spike frequency (Fig. 4B). Electrical activity continues also in our model simulations of P/Qtype channel block with slightly increased spike frequency (Fig. 4C). Reduced Ca 2z entry leads to lower peak Ca 2z concentrations in the submembrane space (½Ca 2z mem ; Fig. 4D). As a consequence, less hyperpolarizing SK-current is activated (Fig. 4E), which leads to an increase in spike frequency (Fig. 4C). Hence, the reduction in excitability caused by blockage of the P/Q-type Ca 2z -current can be overruled by the competing increase in excitability due to the smaller SK-current. Experimentally, v-agatoxin IVA application reduced the action potential amplitude slightly in 3 of 4 cells (by 2.0-4.3 mV), a finding that was quantitatively reproduced by the model, although the reduction was larger (,7.5 mV in Fig. 4C). A direct conclusion from Fig. 4B is that the P/Q-type Ca 2z -current is not needed for the action potential upstroke, unlike the L-type current, probably because of the fact that P/Q-type channels activate at higher membrane potentials than L-type channels. The fact that electrical activity persists with P/Q-type Ca 2z -channels blocked, albeit with lower peak ½Ca 2z mem , could underlie the finding that v-agatoxin IVA only partly inhibits insulin secretion [3].

Paracrine effects on electrical activity
The neurotransmitter c-aminobutyric acid (GABA) is secreted from pancreatic b-cells, and has been shown to stimulate electrical activity in human b-cells [12]. In human b-cells, GABA activates GABA A receptors, which are ligand-gated Cl { channels, thus creating an additional current. Notably, the Cl { reversal potential in human b-cells is less negative than in many neurons, and positive compared to the b-cell resting potential, which means that Cl { currents, such as the GABA A receptor current, stimulate action potential firing in b-cells. Hence, GABA is a excitatory transmitter in b-cells, in contrast to its usual inhibitory role in neurons. We simulate the addition of GABA by raising the GABA A receptor conductance. In a silent model cell with a rather large K(ATP)conductance, simulated GABA application leads to a single action potential whereafter the membrane potential settles at ,245 mV (Fig. 5A), in close correspondence with the experimental results [12]. In an active cell, the simulation of activation of GABA A receptors leads to a minor depolarization and increased action potential firing (Fig. 5B), as found experimentally [12].
Another neurotransmitter, acetylcholine, might also play a paracrine role in human pancreatic islets, where it is released from a-cells, and activates muscarinic receptors in b-cells [22]. Muscarinic receptor activation by acetylcholine triggers a  voltage-insensitive Na z -current in mouse pancreatic beta-cells [23], and similarly, the muscarinic agonist carbachol activates nonselective Na z leak channels (NALCN) in the MIN6 b-cell line [24]. Based on these findings, it was speculated that muscarinic activation of NACLN currents in human b-cells might participate in the positive effect of acetylcholine and carbachol on insulin secretion [4]. Experimentally, we found that carbachol (20 mM) accelerates action potential firing (Fig. 5C). We tested the hypothesis of a central role of leak current activation by increasing the leak conductance in the model to simulate carbachol application, which caused accelerated action potential firing. The simulation thus reproduced the experimental data, and lends support to the hypothesis that carbachol and acetylcholine can accelerate action potential firing via muscarinic receptor-dependent stimulation of NALCN currents [4].

Slow oscillations
We finally use our model to address the origin of slow rhythmic patterns of electrical activity in human b-cells (Fig. 6A) [4,25], which likely underlie slow oscillations in intracellular Ca 2z [26,27] and pulsatile insulin release [28,29]. Based on accumulating evidence obtained in rodent islets [5,30], we have previously speculated that oscillations in metabolism could drive these patterns [7]. In support of this hypothesis, oscillations in ATP levels with a period of 3-5 minutes have been observed in human b-cells [13,14]. By adding a glycolytic component [11], which can oscillate due to positive feedback on the central enzyme phosphofructokinase (PFK), our model can indeed simulate such periodic modulation of the electrical pattern, where action potential firing is interrupted by long silent, hyperpolarized periods, which drives slow Ca 2z oscillations (Fig. 6).

Discussion
Human b-cells show complex and heterogeneous electrophysiological responses to ion channel antagonists. It can therefore sometimes be difficult to reach clear conclusions regarding the participation of certain ion channels in the various phases of parameters except g CaL~0 :100 nS/pF. With default K(ATP)-channel conductance g KATP~0 :010 nS/pF, the simulation reproduces the data in panel C (black traces). When g KATP~0 :002 nS/pF, the model shows continued firing with Na z -channel block (gray traces). E: TTX changed spiking into rapid bursting electrical activity in this human b-cell. F: Simulation showing V (upper), ½Ca 2z mem (middle), and I BK (lower), reproducing the data in panel E. Parameters took default values, except g Na~0 :7 nS/pF, t hNa~3 ms, g Kv~0 :25 nS/pF, g SK~0 :023 nS/pF, g leak~0 :012 nS/pF, and nx PQ~{ 10 mV. The extracellular glucose concentration was 6 mM in all experiments. Each couple of experimental traces (panels A, C and E) is from the same human b-cell before (left) and during (right) application of TTX. In the simulations, the Na z -channel conductance g Na was set to 0 nS/pF during the period indicated by the gray bars. doi:10.1371/journal.pcbi.1003389.g003 We have here shown how mathematical modeling can help in interpreting various electrophysiological responses, and in particular, to study the effect of competing effects and cell heterogeneity. The role of SK-channels in human b-cells is still not clear. We (Fig. 1) and others [8] have found heterogeneous electrophysiological responses to SK-channel antagonists. Our model suggests that these differences can be caused by underlying variations in cell excitability: Less excitable b-cells that produce action potentials evoked mostly by stochastic channel dynamics show a clear increase in action potential frequency when SK-channels are blocked (Fig. 1DE). In contrast, spiking electrical activity in very active cells is driven by the deterministic dynamics caused by ion channel interactions, and is nearly unchanged by SK-channel blockers (Fig. 1A-C). We showed also that rapid bursting activity can be driven by Ca 2z and SK-channels (Fig. 2), which could add a complementary mechanism to HERG-channel dynamics [7] for the control of rapid bursting.
The wide range of responses to TTX could be accounted for by a single model but with different parameters, i.e., differences in the relative size of the various currents. A peculiar finding is the qualitative change from spiking to rapid bursting seen in some cells (Fig. 3E). We suggest that this happens in human b-cells with large Na z -currents. The blockage of this depolarizing current reduces the amplitude of the action potentials, and as a consequence, the size of the hyperpolarizing BK-current. Under the right conditions, the combination of these competing events allows the membrane potential to enter a bursting regime controlled by SKand/or HERG-channels (Fig. 3F). Interestingly, it has been found that TTX reduces insulin secretion evoked by 6 mM glucose greatly, but at glucose levels of 10-20 mM, the effect of TTX on secretion is smaller [3,9,10]. Based on our simulations showing that less excitable cells cease to fire in the presence of TTX (Fig. 3D, black traces), but that lower g KATP can reintroduce spiking activity (Fig. 3D, gray traces), we suggest that at low, nearthreshold glucose levels TTX abolishes electrical activity in many cells, which reduces the ½Ca 2z mem and consequently insulin secretion greatly (Fig. 3D, black traces). At higher glucose concentrations, b-cells have lower K(ATP)-conductance and in some of the cells that stop firing in low glucose concentration the effect of TTX on electrical activity and ½Ca 2z mem is smaller (Fig. 3D, gray traces). Hence, more b-cells remain active in the presence of TTX at high than at low glucose levels. Consequently, insulin secretion is more robust to TTX at higher glucose concentrations.
Similarly, insulin release is more affected by the P/Q-type Ca 2z -channel blocker v-agatoxin IVA at 6 mM (271%) than at 20 mM (231%) glucose [3]. This is in contrast to L-type Ca 2zchannel antagonists, which abolish insulin secretion at both high (15-20 mM) and low (6 mM) glucose concentrations [1,3,10].  Fig. 7A in [12]). Default parameters except g KATP~0 :021 nS/pF. GABA application was simulated by setting g GABAR to 0.1 nS/pF during the period indicated by the gray bar. B: Simulation of application of 10 mM GABA to an active cell (reproducing Fig. 7B in [12]). GABA application was simulated by setting g GABAR to 0.020 nS/pF during the period indicated by the gray bar. Other parameters took default values. C: Experimental recording of spiking electrical activity in the same human b-cell before (left) and during (right) application of carbachol (20 mM). D: Simulation of accelerated action potential firing due to carbachol application. Default parameters except g KATP~0 :016 nS/pF. Carbachol application was simulated by increasing g leak to 0.030 nS/pF during the period indicated by the gray bar. doi:10.1371/journal.pcbi.1003389.g005 These results concerning L-type Ca 2z -channel block are easily explained by the fact that L-type channel activity is necessary for action potential generation [3] (Fig. 4A). In contrast, we showed that electrical activity in human b-cells not only persists, but is accelerated by v-agatoxin IVA (Fig. 4B). The counter-intuitive finding of increased excitability and electrical activity when the depolarizing P/Q-type Ca 2z -current is blocked by v-agatoxin IVA can be accounted for by an even greater reduction in the hyperpolarizing SK-current due to reduced Ca 2z -influx and consequently lower ½Ca 2z mem .
Our mathematical modeling confirmed that GABA released from human b-cells can have a role as a positive feedback messenger. GABA application has been shown to depolarize both silent and active human b-cells [12], which was reproduced here. A detailed characterization of GABA A receptor currents would refine the analysis presented here.
Data from mouse b-cells [23] and the MIN-6 b-cell line [24] suggest that muscaric agonists such as carbachol and acetylcholine stimulate insulin secretion partly by activating NACLN currents. Using our model we could translate this finding to the human scenario, thus testing the hypothesis that this mechanism is also operating in human b-cells [4]. Our simulations confirmed that increased leak currents can underlie the change in electrical activity found experimentally (Fig. 5C). The incretin hormone glucagon-like peptide 1 (GLP-1) has also been shown to act partly via activation of leak channels [31], a mechanism which might be involved in activating otherwise silent b-cells [7,32,33]. These results suggest that leak currents could play important roles in controlling electrical activity in b-cells, and potentially be pharmacological targets. Further studies are clearly needed to investigate these questions.
We were also able to simulate slow rhythmic electrical activity patterns by adding an oscillatory glycolytic component to the model. To date, there is to our knowledge no evidence of oscillations in glycolytic variables in human (or rodent) b-cells or islets, but ATP levels have been found to fluctuate rhythmically also in human b-cells [13,14], supporting the idea of metabolism having a pacemaker role. In agreement, data from rodent b-cells show accumulating evidence for oscillations in metabolism playing an important role in controlling pulsatile insulin secretion [5,30]. It will be interesting to see if these findings in rodents are applicable to human b-cells.
Regarding the model development, the inclusion of SKchannels in the model provided insight that was not within reach with the previous version of the model [7]. Besides the direct investigation of the role of SK-channels, the acceleration in action potential firing seen with P/Q-type Ca 2z channel blockers (Fig. 4) can not be reproduced by the older version of the model without SK-channels [7]. Moreover, considering the effect of TTX on spike amplitude, a better correspondence between experiments and simulations was found with SK-channels included in the model. To model SK-channel activation accurately, we made a special effort to describe ½Ca 2z mem carefully. Submembrane Ca 2z responds rapidly to an action potential, while ½Ca 2z c integrates many action potentials. The rapid submembrane dynamics has important consequences for the study of the role of SK-channels in spiking electrical activity, e.g., it was crucial for explaining the larger effect of TTX on spike amplitude in this version of the model. Most models of electrical activity in rodent b- Figure 6. Metabolically driven slow waves of electrical activity and Ca 2z oscillations. A: Experimental recording of slow oscillations in action potential firing in a human b-cell exposed to 10 mM glucose. B-D: Simulation of slow bursting driven by glycolytic oscillations with glucose concentration G~10 mM and default parameters, except g Kv~0 :2 nS/pF, g SK~0 :02 nS/pF, g BK~0 :01 nS/pF. Oscillations in glycolysis create pulses of FBP (B), which via ATP production modulates K(ATP) channels in a periodic fashion (C). The rhythmic changes in K(ATP) conductance drives slow patterns of electrical activity (D), which causes oscillations in the intracellular Ca 2z concentration (E). doi:10.1371/journal.pcbi.1003389.g006 cells do not include a submembrane Ca 2z compartment, but these models were typically built to explain the slow bursting patterns seen in rodent islets with a period of tens of seconds. For these long time scales, the rapid dynamics in the submembrane compartment is not important. In contrast, the situation is different in human bcells with their faster dynamics.

Modeling
We build on the previously published Hodgkin-Huxley type model for human b-cells [7], which was mainly based on the results of Braun et al. [3], who carefully assured that investigated human islet cells were b-cells. We include SK-channels in the model. Since these channels are Ca 2z -sensitive and located at some distance from Ca 2z -channels [34] we also model Ca 2zdynamics in a submembrane layer controlling SK-channel activity.
The membrane potential V (measured in mV) develops in time (measured in ms) according to dV dt~{ (I SK zI BK zI Kv zI HERG zI Na zI CaL zI CaPQ zI CaT zI KATP zI leak zI GABAR ): ð1Þ All currents (measured in pA/pF), except the SK-current I SK and the GABA A receptor mediated current I GABAR , are modeled as in [7]. Expressions and parameters are given below. For the stochastic simulation in Fig. 1E, we included ''conductance noise'' [35] in the K(ATP) current by multiplying I KATP by a stochastic factor (1z0:2C t ), where C t is a standard Gaussian white-noise process with zero mean and mean square SC t ,C s T~d(t{s), see also [36][37][38]. SK-channels are assumed to activate instantaneously in response to Ca 2z elevations at the plasma membrane but away from Ca 2z channels [34], and are modeled as [39] I SK~gSK Ca n m K n SK zCa m n (V {V K ): ð2Þ In human b-cells, flash-released Ca 2z triggered a ,10 pA current at a holding current of {60 mV, presumably through SK-channels [8].
Assuming that SK-channels were nearly saturated by Ca 2z , the maximal SK-conductance is estimated to be g SK &10 pA= ({60 mV{V K )=C m &0:1 nS/pF. Here, C m = 10 pF is the capacitance of the plasma membrane [3]. In Eq. 2, Ca m is the submembrane Ca 2z concentration (½Ca 2z mem ; measured in mM), which is described by a single compartment model [21] dCa m dt~f where f~0:01 is the ratio of free-to-total Ca 2z , a~5:18|10 {15 mmol/pA/ms changes current to flux, and Vol m and Vol c are the volumes of the submembrane compartment and the bulk cytosol, respectively. B describes the flux of Ca 2z from the submembrane compartment to the bulk cytosol, J PMCA is the flux through plasma membrane Ca 2z -ATPases, and J NCX represents Ca 2z flux through the Na z -Ca 2z exchanger. Cytosolic Ca 2z (Ca c ; measured in mM) follows where J SERCA describes SERCA pump-dependent sequestration of Ca 2z into the endoplasmic reticulum (ER), and J leak is a leak flux from the ER to the cytosol. Expressions and parameters for the Ca 2z fluxes are taken from [40]. The submembrane compartment volume is estimated based on the considerations of Klingauf and Neher [41], who found that a shell model (in contrast to a domain model) describes submembrane Ca 2z satisfactorily when the shell-depth is chosen correctly. The Ca 2z dynamics between channels can be estimated from a shell model at a depth of ,23% of the distance to a Ca 2zchannel. In mouse b-cells the interchannel distance has been estimated to be *1200 nm [42]. Moreover, SK-channels are located w50 nm from Ca 2z channels [34].
Based on these considerations, we modeled the submembrane space controlling SK-channels as a shell of depth *190 nm. The radius of a human b-cell is *13 mm, which gives cell volume (Vol c ), shell volume (Vol m ) and internal surface area (A m ) of the shell, of Vol c~1 :15 pL~1150 mm 3 , Vol m~0 :1 pL, A m~5 30 mm 2 : ð5Þ The flux-constant B can then be calculated as [43] where d m is a typical length scale. We set d m to 1 mm, which together with the diffusion constant for Ca 2z , D Ca~2 20 mm 2 s [41,44], gives B~0:1 ms {1 .
In human b-cells, GABA activates GABA A receptors, which are ligand-gated Cl { channels. We model the current carried by GABA A receptor as a passive current with the expression where g GABAR is the GABA A receptor conductance, and V Cl~{ 40 mV is the chloride reversal potential [4]. We estimate g GABAR from the findings that 1 mM GABA evokes a current of 9:4 pA/pF (but with substantial cell-to-cell variation) at a holding potential of 270 mV [12], which yields a conductance of *0:3 nS/pF. To simulate the changes in firing patterns evoked by lower GABA concentrations (10 or 100 mM) [12], we take into consideration the does-response curve [45] for the a2 b3 c2 subunits, which are the most highly expressed subunits in human b-cells [12]. At 10 mM the GABA-evoked current is w10-fold smaller compared to 1 mM GABA, and we set g GABAR~0 :02 nS/ pF. At 100 mM, the reduction is about 2-fold compared to 1 mM. We used g GABAR~0 :10 nS/pF to simulate application of 100 mM GABA.
To investigate slow electrical patterns (Fig. 6) we added a glycolytic component [11], which drives ATP levels and K(ATP) channel activity. The glycolytic subsystem can oscillate due to positive feedback on the enzyme phosphofructokinase (PFK) from its product fructose-1,6-bisphosphate (FBP). The glycolytic equations are where V GK is the rate of glucokinase, which phosphorylates glucose to glucose-6-phosphate (G6P). G6P is assumed to be in equilibrium with fructose-6-phosphate (F6P), the substrate for PFK, and G6P : F 6P is the sum of G6P and F6P. V PFK is the rate of PFK producing FBP, which is subsequently removed by fructose-bisphosphate aldolase (FBA), which produces glyceraldehyde-3-phosphate (G3P) and dihydroxyacetone-phosphate (DHAP) with rate V FBA . DHAP and G3P are assumed to be in equilibrium, and DHAP : G3P indicates their sum. Finally, G3P serves as substrate for glyceraldehyde-3-phosphate dehydrogenase (GAPDH with rate V GAPDH ), which via the lower part of glycolysis eventually stimulates mitochondrial ATP production. We introduce a phenomenological variable a that mimics ATP levels, and is model by The K(ATP) conductance depends inversely on a, and is modeled as Expressions and parameters are given below. Simulations were done in XPPAUT [46] with the cvode solver, except the stochastic simulation in Fig. 1E, which was performed with the implicit backward Euler method. Computer code can be found as supplementary material, or downloaded from http:// www.dei.unipd.it/ pedersen.

Experiments
Human pancreatic islets were obtained with ethical approval and clinical consent from non-diabetic organ donors. All studies were approved by the Human Research Ethics Board at the University of Alberta. The islets were dispersed into single cells by incubation in Ca 2z free buffer and plated onto 35 mm plastic Petri dishes. The cells were incubated in RPMI 1640 culture medium containing 7.5 mM glucose for .24 h prior to the experiments. Patch-pipettes were pulled from borosilicate glass to a tip resistance of 6-9 MV when filled with intracellular solution. The membrane potential was measured in the perforated-patch whole-cell configuration, using an EPC-10 amplifier and Patchmaster software (HEKA, Lambrecht, Germany). The cells were constantly perifused with heated bath solution during the experiment to maintain a temperature of 31{33 0 C. The extracellular solution consisted of (in mM) 140 NaCl, 3.6 KCl, 0.5 MgSO 4 , 1.5 CaCl 2 , 10 HEPES, 0.5 NaH 2 PO 4 , 5 NaHCO 3 and 6 glucose (pH was adjusted to 7.4 with NaOH). The pipette solution contained (in mM) 76 K 2 SO 4 , 10 KCl, 10 NaCl, 1 MgCl 2 , 5 HEPES (pH 7.35 with KOH) and 0.24 mg/ml amphotericin B. b-cells were identified by immunostaining (18 out of 28 cells) or by size when immunostaining was not possible (cell capacitance .6 pF, [3]). Tetrodotoxin (TTX) and v-agatoxin IVA were purchased from Alomone Labs (Jerusalem, Israel), UCL-1684 was obtained from R&D Systems (Minneapolis, MN), TRAM-34 from Sigma-Aldrich (Oakville, ON, Canada). Figures with experimental responses to ion channel antagonists (Figs. 1, 3, 4 and 5) show recordings from the same cell before (ctrl) and after application of the blocker.

Model equations and parameters
For completeness, we report all expressions and parameters of the mathematical model here. For details, please refer to the Modeling section above and the previous article [7].