To clarify the mechanisms underlying the pancreatic β-cell response to varying glucose concentrations ([G]), electrophysiological findings were integrated into a mathematical cell model. The Ca2+ dynamics of the endoplasmic reticulum (ER) were also improved. The model was validated by demonstrating quiescent potential, burst–interburst electrical events accompanied by Ca2+ transients, and continuous firing of action potentials over [G] ranges of 0–6, 7–18, and >19 mM, respectively. These responses to glucose were completely reversible. The action potential, input impedance, and Ca2+ transients were in good agreement with experimental measurements. The ionic mechanisms underlying the burst–interburst rhythm were investigated by lead potential analysis, which quantified the contributions of individual current components. This analysis demonstrated that slow potential changes during the interburst period were attributable to modifications of ion channels or transporters by intracellular ions and/or metabolites to different degrees depending on [G]. The predominant role of adenosine triphosphate–sensitive K+ current in switching on and off the repetitive firing of action potentials at 8 mM [G] was taken over at a higher [G] by Ca2+- or Na+-dependent currents, which were generated by the plasma membrane Ca2+ pump, Na+/K+ pump, Na+/Ca2+ exchanger, and TRPM channel. Accumulation and release of Ca2+ by the ER also had a strong influence on the slow electrical rhythm. We conclude that the present mathematical model is useful for quantifying the role of individual functional components in the whole cell responses based on experimental findings.
The pancreatic β cell has a unique function of converting variations in the extracellular glucose concentration ([G]) to electrical activity, thereby controlling the level of insulin secretion. This signal transduction is dependent on the interaction between energy metabolism and membrane excitation. Several mechanisms have been suggested underlying this bilateral coupling in pancreatic β cells. The gating of ATP-sensitive K+ channels is regulated by fluctuations in the intracellular concentration of ATP or MgADP ([ATP] or [MgADP]), resulting in a prolongation of the duration of the burst of action potentials with increasing [G]. The activation of L-type Ca2+ channels by an increase of [ATP] (Smith et al., 1989), or the depression of Na+/K+ pump (NaK) activity up to 50% by increasing [G] (Owada et al., 1999), may also favor burst prolongation. In addition, variations in intracellular ion concentrations may have varying influences on individual channels or transporters depending on [G]. For example, it has been recently suggested that a K+ current activated by intracellular Ca2+ (IKslow) may affect bursting activity (Göpel et al., 1999a; Goforth et al., 2002). Finally, the electrical activity induces a significant increase in ion fluxes across the surface membrane, which alters energy consumption via active ion transport or Ca2+-mediated processes, including insulin secretion. These pathways are all linked in a complex system, and one approach to aid the quantification of the contribution to bursting activity of individual pathways is the development of a mathematical β-cell model.
Such models have been used for nearly 30 years to elucidate the principle mechanisms underlying the bursting activity in β cells. Early stage models used a formulation consisting of a minimum number of components: two or three K+ currents, a Ca2+ current, and/or a leak current (Chay and Keizer, 1983; Sherman et al., 1988, 1990; Keizer and Magnus, 1989; Smolen and Keizer, 1992; Bertram et al., 1995b). These model simulations suggested consistently the critical role of a slowly changing variable in generating the burst–interburst rhythm. Subsequent models elaborated metabolic components by including details of glycolysis, tricarboxylic acid (TCA) cycle, and oxidative phosphorylation (Magnus and Keizer, 1998; Bertram et al., 2004; Diederichs, 2006) to examine the gating of IKATP by time-dependent changes in [ADP] or glycolytic oscillation. Several models with detailed descriptions of many more membrane currents and associated changes in intracellular ion concentrations have also been published (Miwa and Imai, 1999; Fridlyand et al., 2003; Meyer-Hermann, 2007).
The object of this study is to clarify quantitatively the detailed ionic mechanisms underlying glucose-induced electrical bursting activity observed in isolated β cells. To achieve this aim, we have developed a comprehensive model based on recent extensive experimental findings on ion channels, transporters, and intracellular Ca2+ dynamics in β cells. If adequate mathematical analyses are successfully applied to this detailed model, the role of individual ion channels will be clarified in quantitative terms, in relation to the principle mechanisms deduced from the theoretical studies using simplified models, and also in relation to the detailed experimental studies on the role of individual functional molecules in real cells.
MATERIALS AND METHODS
The present model of a single β cell was constructed on the framework developed by Fridlyand, Philipson, and their colleagues, the FP model (Fridlyand et al., 2003, 2005), which was designed to examine interactions among glucose metabolism, Ca2+ dynamics including ER, and membrane excitation. The metabolic elements of the model were adopted after minor modifications, whereas the formulations of individual ion channels and transporters were largely revised to reproduce the detailed characteristics of electrical activities reported in the literature. The structure of the model is illustrated in Fig. 1. Because electrical activities or glucose sensitivities vary diversely among different studies or species, we have concentrated on data obtained from dissociated mouse β cells at physiological temperature (33–37°C). Experimental results from other species, including rat and human, or obtained at room temperature were also referred in the absence of relevant mouse data. All equations and parameters are presented in the supplemental material.
Cell dimensions and Ca2+ buffer
Cytosolic (764 fl) and ER (280 fl) volumes and membrane capacitance (6.158 pF) were defined as in the FP model (Fridlyand et al., 2003). This capacitance is within the experimental range measured in isolated β cells in mouse (5.4 ± 0.9 pF) (Rorsman and Trube, 1986) and similar to that found in humans (6.2 ± 0.8 or 7.3 ± 0.4 pF) (Kelly et al., 1991). The concentrations of free Ca2+ in the cytosol ([Ca2+]i) and ER ([Ca2+]ER) were calculated using the buffering power coefficients fi and fER, respectively (see Eqs. S5 and S6 in the supplemental material). The value of fER was determined using the decay time course of the Ca2+ transient evoked by a 45-mM K+ pulse (Gilon et al., 1999), assuming an ER Ca2+-binding capacity of 98–99% as found in gonadotropes (Tse et al., 1994).
Modeling ion channels and transporters
Plasma membrane ion transport comprised eight ion channels and three ion transporters as indicated in Fig. 1. Modeling parameters were based on voltage-clamp data obtained mainly in dissociated β cells, as indicated.
Voltage-dependent Ca2+ current (ICaV).
It is well established that the maximum rate of rise of action potential is determined by the activation of voltage (V)-dependent Ca2+ currents in pancreatic β cells (Ribalet and Beigelman, 1980; Rorsman and Trube, 1986; Ashcroft and Rorsman, 1989). Several types of voltage-gated Ca2+ channels (L, R, and possibly P/Q type) have been reported in mouse β cells (Schulla et al., 2003). Braun et al. (2008) also demonstrated that human β cells express L, T, and P/Q types, but not R-type Ca2+ channels. At present, it is not possible to describe quantitatively each component because of the lack of detailed voltage-clamp data from isolated cells. However, the whole cell Ca2+ currents so far reported show common characteristics for L-type Ca2+ current established in other cell types. These include Ca2+-mediated inactivation (Plant, 1988; Satin and Cook, 1989; Kelly et al., 1991), ultraslow V-dependent inactivation (Satin and Cook, 1989; Kelly et al., 1991), and activation by intracellular ATP (Smith et al., 1989) or “washout” with an ATP-free pipette solution (Hiriart and Matteson, 1988). All of these properties might be heavily involved in modulation of membrane excitability. For example, Henquin and Meissner (1984a) ascribed a gradual decrease in the amplitude and frequency of Ca2+ spikes to Ca2+- and ultraslow V-dependent inactivation of ICaV during the burst. As an initial approximation, we have used a lumped Ca2+ current with characteristics similar to those of L-type Ca2+ current described by a formulation of ICaL developed in cardiac myocytes (Takeuchi et al., 2006). The parameters were adjusted according to voltage-clamp experiments in an insulin-secreting cell line (Satin and Cook, 1989; see Fig. S1) and in isolated mouse β cells (Houamed et al., 2010). Because the kinetics of Ca2+ current of human β cells are similar to those in rodent (Kelly et al., 1991), human data were also used for determining the voltage-dependent activation curve and the time course of ultraslow inactivation. Ca2+-dependent inactivation was described as a function of the single-channel current (Eqs. S18 and S19). For ultraslow inactivation, a noninactivating fraction of 0.4 (Eq. S13) was assumed to reproduce the inactivation curve obtained with a 10-s conditioning pulse (Fig. S1 F). Ideally, all of these properties might better be systematically analyzed in a single study at the physiological temperature. Such a study is awaited.
Delayed rectifier K+ current (IKDr).
Pancreatic β cells show a marked delayed outward K+ current on depolarization (Cook and Hales, 1984; Rorsman and Trube, 1986; Smith et al., 1990b). The voltage-clamp record of IKDr in dissociated β cells of mouse (Houamed et al., 2010) was reproduced (Fig. S2, A and B). For a more precise description of the activation gate (Eqs. S26 and S27), the I-V relationship and activation curve obtained from human (Kelly et al., 1991) and mouse β cells (Rorsman and Trube, 1986) were also considered (Fig. S2 C). In addition, a slow inactivation of IKDr has been observed in various studies (Rorsman and Trube, 1986; Kelly et al., 1991; Houamed et al., 2010). Because the inactivation kinetics were highly dependent on temperature (MacDonald et al., 2003), a time constant of ∼0.3 s at +10 mV at 32–35°C (Houamed et al., 2010) was used for model adjustment (Eqs. S29 and S30).
Ca2+-activated nonselective cation current (ITRPM).
Sturgess et al. (1987) recorded single-channel currents of a Ca2+-activated nonselective cation channel in the INS-1 cell line. Recently, an analogous but more specific current has been described as a TRPM4 channel current, which might be involved in insulin secretion (Cheng et al., 2007; Marigo et al., 2009). In spite of the recent findings, this channel has not been implemented in previous β-cell models. We described the activation by Ca2+ (Eq. S48) with a half-saturation concentration (0.76 µM) and a Hill coefficient (1.7), consistent with experimental recordings (Marigo et al., 2009). The relative permeabilities of Na+ and K+ (Eqs. S50 and S51) were adjusted to give a reversal potential of ∼0 mV (Colsoul et al., 2010). The whole cell conductance of ITRPM was adjusted to reconstruct the plateau potential of approximately −50 mV during the burst.
Store-operated current (ISOC).
Worley et al. (1994a,b) showed that depletion of ER Ca2+ store by zero Ca2+ bath solution with EGTA depolarized the membrane in freshly isolated mouse β cells. They also demonstrated that a nonselective cation current was activated by maitotoxin. Leech and Habener (1998) also recorded a similar maitotoxin-sensitive current that showed a reversal potential of −1.7 mV in insulinoma cell lines. The authors suggested that this current might play a critical role in setting the membrane potential (Vm) to be less negative than the K+ equilibrium potential. Unfortunately, it is still unknown if the store-operated cation currents are attributable to a single class of ion channels at the molecular level. Moreover, the critical level of ER Ca2+ depletion for the half-activation of the store-operated currents (K0.5,ER) remains unclear in β cells. Although a Ca2+ release–activated nonselective cation (CRAN) current has been implemented in previous β-cell models (Bertram et al., 1995a; Chay, 1996, 1997; Mears et al., 1997; Fridlyand et al., 2003), different values of K0.5,ER ranging from 3 to 200 µM were used. In our model, the minimum level of 3 µM [Ca2+]ER is tentatively assumed for the half-activation (Eq. S44). With this assumption, the amplitude of ISOC is minimum under the physiological conditions, but it is activated when the ER is almost completely depleted, for example, by applying thapsigargin.
Miura et al. (1997) demonstrated that depletion of Ca2+ store by thapsigargin triggered Ca2+ influx independent of ICaV, which might be attributed to a different type of current from ICRAN. Based on their finding, ISOC in our model is partly carried by Ca2+, and the size of the current was determined by the steady-state level of [Ca2+]i under thapsigargin and D600, a blocker of ICaV (Eq. S47). It is consistent with the fact that ubiquitous Ca2+ release–activated Ca2+ (CRAC) channels are selective to Ca2+ under physiological ionic conditions (Hoth and Penner, 1993; Prakriya and Lewis, 2002). This Ca2+ entry prevents ER from a serious depletion at 0 mM [G] in our model, and the Na+ and K+ conductance of ISOC sets the resting membrane potential at approximately −70 mV in competition with the background of IKATP (Eqs. S45 and S46).
V- and [Ca2+]i-dependent transient outward K+ current (IKCa(BK)).
A large-conductance Ca2+-activated K+ (BK) current has been recorded in single-channel recordings in insulin-secreting cell lines and mouse β cells (Velasco and Petersen, 1987; Satin et al., 1989; Kukuljan et al., 1991; Houamed et al., 2010). Smith et al. (1990b) found that the amplitude of the whole cell outward current was not affected by chelation of intracellular Ca2+ by adding [EGTA] to pipette solutions, but they observed that a transient component was depressed by blocking ICaV. Furthermore, single-channel recordings demonstrated that this current was activated immediately after the onset of a depolarizing pulse. These findings suggested that the channel might be functionally coupled to Ca2+ channels rather than to bulk cytosolic [Ca2+]. Similar transient outward currents coupled with ICaV have also been reported in human β cells (Herrington et al., 2005; Braun et al., 2008). Because it is difficult to estimate [Ca2+] near the BK channel molecule, this current was tentatively represented as a V-dependent transient K+ current based on the above properties (Eqs. S32–S37). The rate constants for activation and inactivation were determined based on the measurement in dissociated mouse β cells at 33.5°C (Houamed et al., 2010). It has been suggested that this current is a major determinant of the action potential amplitude (Henquin, 1990; Braun et al., 2008; Houamed et al., 2010; Jacobson et al., 2010). Thus, the conductance of IKCa(BK) was determined to set an action potential peak from −10 to 0 mV (Eq. S31).
ATP-sensitive K+ current (IKATP).
It has been well established that the open probability of ATP-sensitive K+ channel changes depending on the intracellular energy status (Cook and Hales, 1984; Rorsman and Trube, 1985), and thereby IKATP modulates membrane excitability and subsequent insulin secretion in β cells (Larsson et al., 1996). Hopkins et al. (1992) suggested that the channel activity is dependent on ADP level over the concentration range of 10–100 µM, rather than on the ATP/ADP ratio (Dunne and Petersen, 1986; Misler et al., 1986). Based on a reaction scheme with two ADP-binding sites (Hopkins et al., 1992), Magnus and Keizer (1998) proposed a detailed model of IKATP. We adopted this model after a minor modification of dissociation constants for ATP and MgADP according to experimental data (Ashcroft and Kakei, 1989; Hopkins et al., 1992).
Rorsman and Trube (1985) found that the input conductance was ∼0.05 nS (20 GΩ in the input resistance) at 10 mM [G], but it increased to 1.9 ± 0.1 nS/pF when ATP was omitted from the intracellular solution. Subsequently, Smith et al. (1990a) observed a similar increase of input conductance to 5.1 ± 0.9 nS under 0 mM [G], which was almost completely inhibited by tolbutamide, a selective KATP channel blocker. Because these data could constrain the maximum conductance of IKATP (GKATP; Eq. S53) in our model, we simulated a corresponding measurement with a voltage-clamp step from −70 to −80 mV. The whole cell input conductance ranged from 0.048 to 0.068 nS during bursting rhythm at 10 mM [G], which was in good agreement with the experimental value (Rorsman and Trube, 1985). Under zero intracellular ATP, however, the input conductance only increased to 0.9 nS, >90% of which was attributed to KATP conductance. This value was much smaller than the experimental measurements. However, we failed to improve GKATP of the original model of IKATP (Magnus and Keizer, 1998), and this problem was left for future work.
Ca2+-activated K+ current (IKCa(SK)).
In islet preparations, Göpel et al. (1999a) recorded a novel K+ current component (IKslow), which was activated with a slow time constant of ∼2.3 s during a train of depolarizing pulses and deactivated with a time constant of 6.5 s after the pulses. An analogous current was also recorded in dispersed mouse β cells in several studies (Göpel et al., 1999a; Goforth et al., 2002; Zhang et al., 2005; Düfer et al., 2009). The pharmacological and gene knockout studies have suggested that small-conductance KCa (SK) channels might contribute substantially to IKslow (Zhang et al., 2005; Düfer et al., 2009). Supporting this view, isoforms of SK −1 to −4 were found to be expressed at the level of mRNA and protein in mouse β cells (Tamarina et al., 2003; Düfer et al., 2009). Interestingly, Kanno et al. (2002) ascribed ∼50% of the experimental IKslow to IKATP. Thus, we implemented the SK channel current as IKCa(SK) in our new model separately from IKATP. The Ca2+ dependency for activation of IKCa(SK) was adopted from Hirschberg et al. (1998) (Eq. S38). It seems that the activation by Ca2+ of SK current is almost instantaneous, but slow changes in [Ca2+]i and/or the contaminated IKATP component might result in the slow time course of IKslow in experimental recordings.
Background nonselective cation current (IbNSC).
Henquin and Meissner (1984a) showed that the resting membrane potential of β cells is less negative than the K+ equilibrium potential. They attributed this depolarizing effect to a basal membrane Na+ conductance (see also Ashcroft and Rorsman, 1989). It is now well established that this background Na+ conductance includes several types of currents. Nevertheless, a background cation current is still required to establish the resting potential, especially when ICRAN is largely inactivated. Thus, we added such a current, IbNSC, of an unspecified nature. Note that many previous β-cell models also included a background current component (Chay and Keizer, 1983; Chay, 1996; Magnus and Keizer, 1998; Meyer-Hermann, 2007; Fridlyand et al., 2009). IbNSC in this model is permeable to Na+ and K+ with a reversal potential at approximately −20 mV (Eqs. S40–S42). The conductance was adjusted to give both the resting membrane potential and input impedance consistent with experimental measurements at a low [G] (Rorsman et al., 1986; Rorsman and Trube, 1986).
Plasma membrane Ca2+ pump (PMCA) and Na+/Ca2+ exchange (NCX) currents (IPMCA, INaCa).
Ca2+ influx through ICaV is balanced with Ca2+ efflux via IPMCA (PMCA1, 2, and 3) and INaCa (NCX1) (Váradi et al., 1995; Herchuelz et al., 2007). PMCA has one Ca2+-binding site and 1:1 Ca2+/ATP stoichiometry (Brini and Carafoli, 2009). PMCA2 has an apparent Hill coefficient of ∼2 (Caride et al., 2001) and the half-maximal concentration of ∼0.1 µM [Ca2+]i in the presence of calmodulin (Enyedi et al., 1991; Elwess et al., 1997). Based on these findings, IPMCA is expressed by a Hill equation (Eq. S95). In addition, it is known that PMCA exchanges one intracellular Ca2+ for one extracellular H+ (Hao et al., 1994), and we assumed that the excess H+ was instantaneously removed by Na+/H+ exchange. Because Na+/H+ exchange was not included in the present model, the resultant Na+ influx by the functional coupling of PMCA and Na+/H+ exchange was directly included in calculating d[Na+]i/dt (Eq. S3).
The description of INaCa was adopted from a cardiac myocyte model (Takeuchi et al., 2006), which describes time-dependent transitions between different functional states of the NCX molecule (Eqs. S75–S94). The slope conductance of INaCa near the reversal potential was 25.5 pS pF−1 at 14 µM [Ca2+]i and 30 mM [Na+]i in the present model, which is about half of the experimental value (53 pS pF−1) (Gall et al., 1999). This difference seems to fall within the range of experimental variations because of the limited intracellular perfusion with pipette solutions through the ruptured patch.
NaK current (INaK).
The INaK model was adopted from Oka et al. (2010), in which the turnover rate was precisely described in terms of Vm, intracellular, and extracellular compositions of Na+ and K+, and the free energy of ATP hydrolysis (ΔGATP) based on thermodynamics (Eqs. S54–S74). Although this model was developed with reference to experimental measurements in cardiac myocytes, we assumed for convenience that the basic characteristics of the pump activity would be common in β cells. In addition, the inhibition of the pump activity by glucose via intracellular signaling (Owada et al., 1999) was implemented (Fglc; Eq. S55). The amplitude factor of INaK (PNaK) was determined to satisfy Na+ homeostasis in both quiescent and bursting activities. Finally, the K+ balance between efflux through K+ channels and the active influx via NaK was calculated, rather than fixing [K+]i as in the original FP model.
Modeling intracellular Ca2+ dynamics
A precise description of ER Ca2+ dynamics is critical for modeling β-cell function. Uptake of Ca2+ into the ER is mediated by ER Ca2+ ATPase (SERCA), and approximately equal amounts of SERCA 2b and 3 are expressed in pancreatic islets (Váradi et al., 1996). The apparent affinity for cytosolic Ca2+ was determined with a half-activation concentration (K1/2) of 0.27 and 1.1 µM, and a Hill coefficient (nH) of 1.7 and 1.8 for SERCA 2b and 3, respectively (Lytton et al., 1992). The SERCA activity in the present study was represented with a Hill equation of K1/2 = 0.5 µM and nH = 2, compromised for the whole cell simulation (Eq. S96). Ca2+ release from ER is a critical determinant for reconstructing the slow decay phase of [Ca2+]i observed after action potential burst. Although an application of IP3 facilitates Ca2+ release (Tengholm et al., 2001), the slow Ca2+ decay during the interburst did not seem to be triggered by IP3-, depolarization-, nor Ca2+-induced Ca2+ release (Gilon et al., 1999). Therefore, ER Ca2+ release (Jrel) was described as a passive flux down a concentration gradient in this study (Eq. S97).
ER volume (volER), maximum velocity of SERCA (PSERCA; Eq. S96), nor the permeability of the Ca2+ release channel (Prel; Eq. S97) has been fully measured to provide definite values of these parameters. Thus, they were adjusted based on the following experimental findings. (a) The physiological level of [Ca2+]i hardly exceeds 0.5 µM during glucose stimulation (Rorsman et al., 1984). (b) The resting [Ca2+]i is 60–100 nM (Rorsman et al., 1992; Chow et al., 1995). (c) Onset and offset time courses of Ca2+ transient were recorded, which were evoked by the action potential burst, a voltage-clamp pulse, or K+-induced depolarization (Gall et al., 1999; Gilon et al., 1999). (d) Direct measurement of [Ca2+]ER using a low affinity Ca2+ fluorescent dye revealed that [Ca2+]ER is maximally increased up to ∼200 µM by Ca2+ uptake through SERCA in the absence of IP3 (Tengholm et al., 2001). It was consistent with [Ca2+]ER of 60–200 µM suggested previously (Tse et al., 1994). (e) At 12 mM [G], Ca2+-stimulated ATPase activity of SERCA was comparable to that of PMCA in β cells (Roe et al., 1994). In the present β-cell model, the ratio of ATP consumption by SERCA and PMCA was approximately 1:1 at 12 mM [G], ranging from 1:3 in a quiescent state at 6 mM [G] to 4:3 during continuous firing at 20 mM [G].
Modeling energy metabolism
Fridlyand et al. (2005) elaborated a set of equations for ATP production through glycolysis and oxidative phosphorylation, and for ATP consumption based on a wide range of biochemical studies. We used their model with a few modifications as follows. First, we changed the glucose dependency of glycolysis (fglc) (Eq. S100) to reproduce the experimental finding that the burst duration is prolonged with increasing [G] in β cells. Our revision might be appropriate because fglc reflects the [G] dependency of all the reaction steps including glycolysis and TCA cycle in our model. Note that the original values in the FP model were determined under the assumption that glucose phosphorylation by glucokinase was the only limiting step in glycolysis. Second, we calculated ATP production via β oxidation of fatty acid (Jβ,ox; Eq. S99), in addition to glycolysis (Jglc; Eq. S98). This modification prevented the system from a metabolic collapse at a low [G] (<2 mM), which actually occurred in the FP model. Third, in the production of reduced metabolic compounds (Re), we took account of the total amount of pyridine nucleotides ([Retot]) by adding a term of ([Retot] − [Re]) in Jglc and Jβ,ox (Eqs. S98 and S99). This term was crucial to avoid an unlimited increase of [Re] at a high [G] (>15 mM), observed in the FP model. Under the assumption that most Re consists of NADH in the mitochondria, [Retot] of 10 mM was used (Cortassa et al., 2003). The consumption of [Re] by oxidative phosphorylation was calculated using a stoichiometry of 2.5 between ATP and NADH, and with a volume ratio (2.5) between the cytosol and mitochondria (Eq. S102).
Lead potential (VL) analysis
To clarify the ionic mechanisms underlying burst–interburst rhythm in our new β-cell model, we applied the VL analysis developed by Cha et al. (2009). The method quantifies the contributions of individual membrane currents to changes in Vm by calculating an equilibrium potential at each moment (VL) using the time-varying conductance (GX), reversal potential (EX), and V-independent transporter current (IY),
Also refer to Eq. S108. VL always moves in advance of Vm, and its time derivative (dVL/dt) drives the automatic change of Vm. The relative contribution (rc) of a current component of interest (i) is defined by a relative change in dVL/dt when the time-dependent change of i is selectively fixed. The total sum of rc for all components equals unity at each time point, and is used to validate the calculations,
This method has been verified in various cardiac cell models (Cha et al., 2009; Himeno et al., 2011). In the present study, the contribution c (mV s−1) was used, instead of rc. c was newly defined by the following equation:
c with a positive sign indicates that the corresponding component contributes to membrane depolarization, and vice versa.
Among the three electrogenic ion transporters, V-independent IPMCA was treated as a current source (Eq. S108). INaK and INaCa were expressed with Eqs. 4 and 5, where GNaK and GNaCa are the slopes of tangential lines fitted to the instantaneous I-V relation at each moment, and Ex_NaK and Ex_NaCa, the intersections of the tangential lines with the x axis:
The contribution of INaK or INaCa in Fig. 5 was a summation of c evaluated by fixing GNaK and Ex_NaK, or GNaCa and Ex_NaCa, respectively. Because GNaK and Ex_NaK are functions of [Na+]i, [K+]i, [ATP] or [MgADP], and Vm, the contribution of each concentration change was also evaluated in the bottom panels of Fig. 5.
Online supplemental material
Equations, parameters, and the definition of symbols of the β-cell model are provided in the supplemental material. Table S1 lists the initial values of the 18 variables in this model. Figs. S1 and S2 show reconstructions of ICaV and IKDr in voltage-clamp experiments, respectively. Fig. S3 shows the effect of thapsigargin on the Ca2+ transients induced by applying high K+ pulses to the model. Fig. S4 is VL diagram of the FP model for comparison to our model (Fig. 5). The supplemental material is available.
Electrical activity and intracellular concentrations of ions and metabolites in pancreatic β cells
Burst of action potentials evoked by various glucose concentrations.
Fig. 2 shows the time-dependent changes in Vm, [ATP], [MgADP], [Na+]i, and [Ca2+]i evoked by different [G] in the new β-cell model. At [G] < 6 mM, the membrane was quiescent, and the concentrations of intracellular ions and metabolites remained at various steady-state levels depending on [G]. The resting potential decreased from −70 mV at 0 mM [G] (not depicted) to −58 mV at 6 mM [G], accompanied by an increase in the input impedance from 4 to 15 GΩ. This input impedance is comparable to experimental measurements of 3–30 GΩ (Rorsman and Trube, 1986), 1–10 GΩ (Rorsman et al., 1986), or 3 GΩ (Smith et al., 1990a). In the simulation, the increase in input impedance largely resulted from the progressive closure of IKATP channels. At 7 mM [G], a typical burst of action potentials appeared. The burst duration was elongated as [G] increased, and finally the burst was transformed to a continuous firing at [G] > 19 mM (Ashcroft et al., 1984; Henquin and Meissner, 1984b). The interburst phase is also elongated at a higher [G] in the present study because more time was required to recover from ion accumulation during the preceding burst period of longer duration. This simulation result is in agreement with the experimental data from mouse islets showing longer burst and interburst periods at a higher [G] (Antunes et al., 2000). Our model, however, failed to reconstruct gradual shortening of the interburst period with [G] (Meissner and Schmelz, 1974).
The action potential in the model is in good agreement with the representative burst activity recorded in a single β cell in the presence of 2.6 mM [Ca2+]o and 10 mM [G] at 31°C (see Fig. 1 B in Smith et al., 1990a). The maximum rate of rise was 2–3 V s−1 in the model, comparable to 3.2 V s−1 (Rorsman and Trube, 1986) or 3.5 V s−1 (Dean et al., 1975). The peak potential was about −4 mV in the model versus −8.3 mV experimentally (Smith et al., 1990a), the plateau potential was about −50 versus −53.7 mV, and the maximum negative potential during the interburst period was about −68 versus −76.4 mV. The maintenance of the plateau potential was mainly attributable to ICaV conductance remaining at the end of the action potentials. It was supported by a simulation showing that the burst was interrupted if ICaV was instantaneously deactivated by applying a brief hyperpolarizing voltage pulse (not depicted). The Ca2+-activated inward currents, ITRPM and INaCa, also contributed to the maintenance of the plateau potential.
Slow fluctuations in [ATP], [MgADP], [Na+]i, and [Ca2+]i during burst–interburst rhythm.
In our model, [ATP] and [MgADP] changed in synchrony with electrical events at [G] > 7 mM (Fig. 2, second row). That is, [MgADP] increased at the expense of ATP during the burst and in turn decreased during the subsequent quiescent period when the cell was relieved from the extra Ca2+-dependent ATP consumption. These typical responses were observed at 8 mM [G]. At 12 or 16 mM [G], however, the ATP consumption was compensated for to a greater extent by increased ATP production. Thus, [MgADP] increased much slower during the burst, and its maximum level at the end of burst was lower in spite of the elongated burst duration. On the other hand, the fluctuation in [Na+]i was enlarged with an increase in burst duration, and finally [Na+]i remained elevated at [G] > 19 mM (Fig. 2). Accumulation of [Na+]i was mostly a result of Na+ influx through NCX, which compensated for the large Ca2+ influx through ICaV. Based on the opposite changes in the fluctuations of [ATP] and [Na+]i by increasing [G], our β-cell model predicted that the activation of INaK by the accumulation of [Na+]i might take over the role of IKATP in terminating the burst at a higher [G].
Fluctuation in [Ca2+]i during the burst–interburst rhythm also has profound effects on the electrical activity. As demonstrated in Fig. 2, [Ca2+]i jumped from a resting level of ∼100 to ∼400 nM at the onset of the burst, and then the plateau level of the oscillation (fast Ca2+ ripple) slowly decreased during the burst, because of the slow inactivation of ICaV. At [G] > 12 mM, a brief oscillation in the plateau level of the Ca2+ ripple preceded the final termination of the burst, which has not been described by experimental studies. We found that this oscillation was sensitive to the amplitude of IKCa(SK) but failed to clarify the underlying mechanisms in the present study. After cessation of the burst, a slow decay phase (or Ca2+ tail) was observed at 12 and 16 mM [G], but hardly at 8 mM [G]. This Ca2+ tail is caused by release of Ca2+ from the ER, which has accumulated during the preceding burst. The increase in the Ca2+ fluctuation at a higher [G] has complex influences on membrane ion channels or transporters, that is, activation of outward-going IPMCA or IKCa(SK), as well as inward-going INaCa or ITRPM. The overall effects of [Ca2+]i will be evaluated mathematically later.
Role of ER Ca2+ dynamics in glucose-induced burst–interburst rhythm
Ca2+ dynamics in the new β-cell model were validated before we analyzed the ionic mechanisms. In control conditions, a regular burst–interburst rhythm and the accompanying Ca2+ transients were generated with a cycle length of ∼40 s at 11 mM [G] (Fig. 3, the left half). At the onset of a burst, most Ca2+ influx through ICaV was instantaneously captured by cytosolic Ca2+-binding proteins (fi in Eq. S5). Then, during the initial 1 s of the burst, the Ca2+ influx was compensated for by the ER (JSERCA-Jrel; 44%), PMCA (26%), and NCX (34%) (Fig. 3, bottom), which was in good agreement with experimental results (Gall et al., 1999). As the burst progressed, Ca2+ gradually accumulated in the ER, and thus the ER Ca2+-buffering capacity became less effective because of an increase in Ca2+ release from the ER. Importantly, 97% of the Ca2+ accumulated during the whole burst was taken up by the ER, and only 3% remained in the cytosol. After cessation of the burst, the accumulated Ca2+ in the ER was slowly released into the cytosol (Fig. 3, bottom), which is a main contributor of the long-lasting Ca2+ tail. This simulation result is in line with experimental responses (Gilon et al., 1999).
For further examination of the relevance of the Ca2+ dynamics in the model, the effects of blocking SERCA by thapsigargin were simulated. In the right half of Fig. 3 indicated by a gray horizontal bar, the activity of SERCA was reduced to 20% of the control. As a result, the Ca2+-buffering capacity of ER decreased, and in the steady state, the amplitude of Ca2+ oscillation was increased by nearly two times. In addition, the Ca2+ tail disappeared from the interburst period and the electrical rhythm became about two times faster through the shortening of both interburst and burst periods. These findings are in good agreement with several experimental recordings (Miura et al., 1997; Gilon et al., 1999; Fridlyand et al., 2003) and previous simulation results (Fridlyand et al., 2003; Bertram and Sherman, 2004). The rate of depolarization during the interburst was accelerated by the activation of inward ISOC as a result of ER depletion. The burst duration was also reduced because the opening of IKATP was accelerated by the enhanced Ca2+-dependent ATP consumption. Increased outward IKCa(SK) or IPMCA by the amplified Ca2+ transient might also help the early termination of the burst, whereas inward INaCa and ITRPM have the opposite effects.
We also simulated Ca2+ transients induced by applying 45 mM of K+ solution (Fig. S3). The Ca2+ tail observed after the high K+ pulse was well reconstructed (Gilon et al., 1999). The simulation predicted that [Ca2+]ER was accumulated up to ∼60 µM via ICaV activated through high K+-induced depolarization (approximately −25 mV). In the presence of thapsigargin, the amplitude of Ca2+ transients was increased with a large initial peak, and the slow Ca2+ tail disappeared. The slow inactivation of ICaV caused the marked decrease in [Ca2+]i during the initial 10 s of the pulse, as well as the temporal depression after washing out the high K+ solution.
Ionic mechanisms underlying the electrical activity of β cells
Current profile during the burst and interburst periods.
The findings in Fig. 2 suggested that the burst rhythm is determined by the balance among current components that are modulated by slow changes in [ATP] and [MgADP], as well as those in [Na+]i and [Ca2+]i. We measured the amplitudes of all these currents, including IKATP, INaK, INaCa, IPMCA, ITRPM, and IKCa(SK), at 8 and 16 mM [G], in addition to V-dependent ICaV and IKDr (Fig. 4). During the burst period, the current levels were measured at the most negative potential between successive action potentials. The plateau potential gradually shifted negative toward the threshold for the full repolarization of the burst termination. At both [G], IKDr was of minimum size because of almost complete deactivation at the end of individual action potentials, and its contribution to changing the plateau potential seemed to be negligible. In contrast, ICaV had the largest amplitude, suggesting that it is the major current maintaining the plateau potential or driving the interburst depolarization to trigger the subsequent action potential burst. IKATP provided a sizable outward current during the interburst at 8 mM [G] but was much decreased at 16 mM [G]. In contrast, outward INaK and IKCa(SK), and inward INaCa and ITRPM, were substantially increased at 16 mM [G] by the accumulation of [Na+]i and [Ca2+]i during the prolonged burst period. The amplitude of ISOC was negligibly small throughout the records in Fig. 4 at both 8 and 16 mM [G] (not depicted). These current profiles, however, only give clues as to the contribution of individual currents underlying the generation of electrical bursting activity. A quantitative understanding of the ionic mechanisms requires further mathematical analysis, such as VL analysis in the next section or bifurcation analysis as described in our companion paper (see Cha et al. in this issue).
VL analysis of interburst ionic mechanisms.
To measure the contribution of each current component to automatic change in Vm, VL analysis was applied to the simulation results (Eqs. 1, 3, and S108). The magnitudes of the contribution (c in mV s−1; see Materials and methods) of individual ion channels and transporters were calculated over the interburst period, as indicated with horizontal gray bars in Fig. 4 (A and B). c was plotted in a cumulative manner at 8 and 16 mM [G] (Fig. 5, middle panels).
At 8 mM [G], V-dependent activation of ICaV (dCaV), albeit a tiny change from 0.03 to 0.05, provided the largest positive contribution during the entire course of slow depolarization (Fig. 5 A). In contrast, the contribution of ultraslow inactivation of ICaV (fus) was trivial. IKATP, an outward current, also provided a positive contribution to the depolarization (c ∼ 0.1–0.2 mV s−1) because its open probability was gradually reduced by both increasing [ATP] and decreasing [MgADP]. In the late phase, the contribution of IKATP became smaller by gradual equilibration of [ATP] and [MgADP]. The positive contribution of inward INaCa (c < 0.1 mV s−1) was mainly attributable to increased turnover rate by the gradual decrease of [Na+]i after cessation of the burst. INaK, IKCa(BK), and ITRPM hindered the slow depolarization, as represented by their negative contributions (less than −0.2 mV s−1).
At 16 mM [G], the ionic mechanisms changed markedly (Fig. 5 B). The contribution of IKATP almost disappeared from the VL diagram, but the contribution of Ca2+-dependent currents (IPMCA, INaCa, and ITRPM) noticeably increased in compared with those at 8 mM [G]. The Vm change showed two phases during the interburst period: early hyperpolarization and late depolarization. During the early phase, the hyperpolarization was mainly attributed to decreases in inward INaCa and ITRPM as a result of the slow decay of [Ca2+]i. The sum of these hyperpolarizing effects was larger than the depolarizing effect caused by the decrease in outward IPMCA. In the late phase, the decay rate of [Ca2+]i slowed down, the contribution of INaCa was reversed by the decrease in [Na+]i, and the negative contribution of ITRPM was also reduced. Furthermore, the decrease in [Na+]i gradually reduced outward INaK and contributed to depolarization. As a consequence, the membrane started to depolarize at the late phase.
Comparison of the VL diagrams in Fig. 5 (A and B) reveals that a metabolic-dependent mechanism (IKATP) at a lower [G] was replaced by an ion-dependent mechanism (IPMCA, INaCa, and ITRPM) at a higher [G] in generating the burst–interburst rhythm. This replacement of mechanism was further exemplified by separating the contribution of INaK into metabolism- and ion-dependent mechanisms (Fig. 5, bottom panels). At 8 mM [G], a negative contribution of INaK was caused by rapid recovery of the ATP/MgADP composition, whereas at 16 mM [G], the metabolic effects almost disappeared, and the decrease in [Na+]i dominated the time course of c of INaK.
VL analysis of repetitive action potentials.
The result of VL analysis is presented in Fig. 6 for two successive action potentials during the burst. The VL (Fig. 6, red line) leads the time-dependent change in Vm (black line) in advance and intersects the Vm curve when dVm/dt (or Itot) equals zero. The VL diagram (Fig. 6, bottom) indicates that the time course of the action potential is largely determined by ICaV. In the rising phase of the spontaneous action potential, the progressive V-dependent activation of ICaV plays the major role; likewise, the V-dependent deactivation of ICaV is mainly responsible for repolarization. The activation of IKCa(BK) partially counteracts ICaV to reduce the maximum rate of rise or decay of the action potential. Surprisingly, the delayed activation of outward IKDr provided a negative contribution only at the beginning of the repolarizing phase, but then reversed its contribution to retard the repolarizing influence of ICaV. This retarding effect of IKDr is a result of V-dependent removal of activation (pKDr). The contributions of the other substrate-dependent currents, IKATP, IPMCA, INaCa, ITRPM, and INaK, are barely visible because the concentrations of ions or metabolites changed minimally over the time span of an action potential.
An extra effect of [G] on the bursting activity through direct inhibition of NaK
Owada et al. (1999) demonstrated that applying glucose to β cells inhibited Na+/K+ ATPase in a dose-dependent and reversible manner via a distinct signal transduction pathway. Because this inhibition was of considerable magnitude (up to 55%), they suggested that the inhibition of INaK might promote insulin secretion at a high [G]. We tested this hypothesis by switching on the inhibitory action of glucose on INaK (Fglc; Eq. S55) after a steady rhythm was established (Fig. 7). Immediately after INaK was reduced by introducing the glucose inhibition (Fig. 7, gray bar), an action potential burst of longer duration was evoked accompanied by a larger Ca2+ transient, in agreement with the experimental observations using an NaK blocker (Bozem and Henquin, 1988). Contrary to the expectation of Owada et al. (1999), the burst interval returned to control at the next burst and remained constant. The amplitude of INaK was almost restored because [Na+]i gradually increased until INaK exactly matched the Na+ influx. The basal level of [Ca2+]i was initially increased by the intervention but slowly recovered over the next 100 s. Similar results were simulated at 12 mM [G]. The simulation suggests that the partial inhibition of INaK by glucose might increase insulin secretion at 8 mM [G], but the effect is only transitory.
By integrating a broad range of electrophysiological findings into a mathematical model, the response of pancreatic β cells to extracellular glucose was well reconstructed, and the underlying mechanisms were elucidated in a comprehensive manner. The new β-cell model showed a series of responses to varying [G], that is, the intermittent burst of action potentials accompanied by Ca2+ transients at [G] > 7 mM, the elongation of the burst duration with increasing [G], and the continuous firing of action potentials at [G] > 19 mM. VL analysis of the model successfully quantified contributions of ion channels and transporters to the slow interburst depolarization. It was concluded that alternating burst and interburst events at the physiological range of [G] is regulated mainly by IKATP channels, which transduce signals from varying [ATP] or [MgADP] to membrane excitability. The novel prediction is that the role of IKATP is taken over by electrogenic ion transporters, such as INaCa, INaK, IPMCA, and a Ca2+-activated ion channel, ITRPM, at a higher [G].
Comparison with the FP model
To our knowledge, the β-cell model developed by Fridlyand et al. (2003, 2005) provided the first description of individual channels and transporters on a plasma membrane at a molecular level. Our model is based on the structure of this FP model to couple membrane excitation with energy metabolism. We revised most of the ionic current components with reference to more extensive electrophysiological findings. In the FP model, a high K+ external solution induces continuous Ca2+ influx through ICaV (about −30 to −50 pA) and eventually causes a metabolic collapse by a rapid depletion of cytosolic ATP. Relevant simulations to experimental findings were obtained when both Ca2+-mediated inactivation and V-dependent ultraslow inactivation were included in the new model of ICaV. Moreover, we added new currents, ITRPM and IKCa(BK), based on recent experimental findings. We found that ITRPM is an important current to maintain the plateau potential around −50 mV during an action potential burst, whereas a full repolarization between action potentials was observed in the FP model. IKCa(BK) is important in the regulation of action potential amplitude.
For self-consistency of the model, we included all ion transports across the cell membrane in calculating both Vm and intracellular ion concentrations, according to charge conservation law (see Cha et al., 2011). (a) We took account of the H+ influx via Ca2+/H+ exchange through PMCA. This H+ flux was assumed to be completely converted to equivalent Na+ flux by a fast Na+/H+ exchange (see Materials and methods). (b) [K+]i was not fixed in our model, but the time-dependent change was calculated by K+ fluxes through NaK and ion channels. These modifications were prerequisite for examining the roles of ion transporters, which are greatly affected by ion concentrations.
The model of INaK in the FP model was largely modified from the original model of Chapman et al. (1983) by omitting state transitions of the carrier protein. This simplification resulted in no saturation of the turnover rate by [Na+]i. Moreover, INaK in the FP model is a V-independent current in the range from −80 to 0 mV. We implemented a new kinetic model of INaK with the state transitions (Oka et al., 2010). It shows properties well established in experimental studies, such as dependencies on Na+ and K+, the free energy of ATP hydrolysis (ΔGATP), and the membrane potential. Thereby, the present study reliably predicted that INaK took a pivotal role in terminating the burst when [Na+]i was accumulated during a long-lasting burst at a high [G] (Fig. 5). For comparison with the present study, we applied VL analysis to the slow interburst depolarization in the FP model (Fig. S4). At a relatively low [G] (8.5 mM), a VL diagram demonstrated that the time-dependent decrease of outward INaK takes the major role in determining the depolarization rate (c ∼ 0.25 mV s−1), whereas the contribution of IKATP was negligibly small (c ∼ 0.025 mV s−1). It is because of a relatively rapid production of ATP in the FP model, resulting in a long-lasting burst even at relatively lower [G], accompanied by [Na+]i oscillation with an amplitude of ∼2 mM.
Mechanisms to generate the bursting activity in modeling studies
Complex patterns of electrical activity in β cells with varying [G] has been one of the interesting targets in the field of mathematical physiology, and several explicit hypotheses have been put forward in various forms of mathematical models. Importantly, however, the fundamental question still remained as to what the slowly varying factor underlying the time course of burst–interburst rhythm in β cells is. Here, we discuss the multiple key membrane components suggested in relation to the slow intracellular factors hypothesized in previous modeling studies, that is, [Ca2+]i, [Ca2+]ER, [ATP] and/or [ADP], and [Na+]i.
[Ca2+] and the Ca2+-activated K+ currents.
One of the major hypotheses assumed a gradual activation of Ca2+-dependent K+ currents during an action potential burst. For example, early Chay-Keizer models adopted the BK channel (Chay and Keizer, 1983; Sherman et al., 1988), and they predicted that an accumulation of intracellular Ca2+ via repetitive spikes increased the outward BK current and terminated the burst. In turn, the burst was resumed when the BK channels were sufficiently deactivated during the interburst period. Distinct from the expectation in their models, however, the progressive accumulation of Ca2+ has not been established experimentally, but rather, a rapid rise of the Ca2+ transient leveled off to the plateau level within the initial several seconds of the burst (Santos et al., 1991; Worley et al., 1994a), or [Ca2+]i slightly rose (Gilon and Henquin, 1992; Zhang et al., 2003) or decayed (Miura et al., 1997; Henquin et al., 2009; Merrins et al., 2010) thereafter. In recent studies, an existence of a different type of Ca2+-dependent K+ channel, IKslow, has been reported (Göpel et al., 1999a; Goforth et al., 2002; Zhang et al., 2005). The important role of IKslow in terminating the burst via [Ca2+] accumulation was suggested by two mathematical models (Goforth et al., 2002; Fridlyand et al., 2009). To simulate the expected slow kinetics of IKslow, Goforth et al. (2002) assumed a localized subspace with a substantial volume (occupying ∼30% of the cytosolic volume), in which the Ca2+ concentration activating IKslow varies in parallel to changes in [Ca2+]ER, rather than the more rapid variation of [Ca2+]i. However, such kinetic changes in IKslow during the bursting activity nor any histological evidence to justify the diffusion barrier have been found experimentally. Alternatively, Fridlyand et al. (2010) assumed IKCas with an extremely slow time constant of 2.3 s for activation to represent IKslow. However, SK channels, one of the candidates contributing to the experimental IKslow, show very fast gating kinetics (Hirschberg et al., 1998).
The slow decay phase of [Ca2+]i after the burst might lead to slow changes in membrane conductances of several Ca2+-activated channels or transporters. Indeed, the present study suggests that [Ca2+]i might play a significant role in driving slow interburst depolarization through IPMCA, INaCa, and ITRPM, as well as IKCa(SK) (Fig. 5). Unfortunately, only a few indirect measurements of these currents have been reported experimentally. As a result, the contribution of ITRPM was only considered in our model, and INaCa and IPMCA were implemented in a few previous models (Fridlyand et al., 2003; Diederichs, 2006; Meyer-Hermann, 2007).
[Ca2+]ER and ISOC.
Repetitive empting and refilling of ER is closely related to bursting rhythm by modulating Ca2+-activated currents. Among them, a store-operated inward current, ISOC (sometimes termed ICRAC or ICRAN), would be a primary candidate to generate the bursting rhythm. Gilon et al. (1999) suggested that ER fills with Ca2+ during the burst, and the gradual deactivation of ISOC may lead to termination of the burst. Conversely, the subsequent emptying of the ER after the burst might then reactivate ISOC to trigger a new burst. This mechanism has been tested in several models (Bertram et al., 1995a; Chay, 1996, 1997; Mears et al., 1997; Fridlyand et al., 2003). However, the half-activation concentration of [Ca2+]ER (K0.5,ER) has not been measured experimentally, and thus the predicted contributions of ISOC are different among studies. For example, in the models of Chay (1996, 1997) using a K0.5,ER of 50 or 70 µM, the gating of ISOC took a central role in determining the burst rhythm. On the other hand, in the other models ISOC contributed little because the channel remained closed as a result of the assumption of a relatively low K0.5,ER (3 µM; Bertram et al., 1995a, and the present model), or was always open because of an assumed high K0.5,ER (200 µM; Fridlyand et al., 2003). Therefore, it is important for K0.5,ER to be determined experimentally to decide the role of ISOC in generating glucose-induced bursting rhythm under normal conditions. Interestingly, rather consistent effects (the prolongation of the spike burst or the acceleration of the bursting rhythm) were reconstructed with these differing models when ISOC was maximally activated by ER depletion under thapsigargin, muscarinic antagonist, or low glucose (Bertram et al., 1995a; Mears et al., 1997; Fridlyand et al., 2003; the present model).
[ATP] and [ADP] and IKATP.
[ATP] and/or [ADP] have been considered as key slow factors, and several β-cell models examined the time-dependent gating of KATP channels. However, quantitative estimation of the contribution of IKATP to burst activity is highly dependent on the formulation of both IKATP and the metabolic components of each model. For example, Magnus and Keizer (1998) developed a detailed IKATP model and concluded it was a major factor, whereas simulations using the FP model and the same IKATP formulation concluded that IKATP was not of major significance. This is because the two models adopted radically different schemes describing the production of [ATP] and [MgADP].
In addition, it should be noted that ATP-consuming transporters, such as PMCA, SERCA, and NaK, should influence the bursting rhythm by modulating their activities according to the intracellular energy level. However, few studies have dealt with this subject, except the present model by incorporating the detailed kinetic model of INaK with ΔGATP dependency.
[Na+]i and INaK.
As demonstrated here and in previous modeling studies (Miwa and Imai, 1999; Fridlyand et al., 2003; Meyer-Hermann, 2007), glucose-induced fluctuations of ICaV results in rhythmical Na+ entry through the action of NCX. Increased [Na+]i will activate INaK and lead to termination of the burst; in turn, a slow decay of [Na+]i leads to a decrease in INaK during the interburst period. Experimentally, Grapengiesser (1996, 1998) observed distinct oscillations of [Na+]i under a partial suppression of NaK in mouse β cells. In support of this idea, these oscillations disappeared after inhibition of ICaV or under a lower glucose, but they were insensitive to a blocker of V-dependent Na+ channels.
[Na+]i also modulates the turnover rate of the NCX exchanger. Thus, a proper Na+ dependency of INaCa is essential for examination of the role of [Na+]i in generating bursting activity. The kinetic scheme of INaCa used in this study was developed in cardiac cells and has been well tested experimentally. In addition, the tetrodotoxin-sensitive Na+ current (INa) might also contribute to intracellular Na+ accumulation. In preliminary studies, we implemented INa based on recordings in pancreatic β cells from rat (Hiriart and Matteson, 1988), mouse (Göpel et al., 1999b; Vignali et al., 2006), and human (Braun et al., 2008). However, because INa was almost completely inactivated at the physiological Vm, the generated flux was trivial. VL analysis also revealed that INa made a very minor contribution to the slow depolarization. Thus, INa was deleted from the present model.
The [Na+]i of 10–14 mM recorded by Grapengiesser (1996) is much higher than the 5.5–7.5 mM in our simulations. When examined with our model, however, 10–14 mM [Na+]i resulted in the reverse mode of the NCX all through the normal burst activity because the Na+-driving force is much reduced. Furthermore, the amplitude of the oscillation of [Na+]i and the corresponding effect on INaK were reliably estimated in our study. This is because the average Na+ influx through the NCX was determined by the amplitude of ICaV and the action potential frequency, both of which were consistent with experimental data.
Is a single β cell capable of generating full-sized action potentials?
Remarkably, the action potential parameters in our model are quite comparable to experimental measurements in isolated mouse β cells obtained by Smith et al. (1990a) (see Slow fluctuations in [ATP]… in Results). In most papers, however, the amplitude of action potentials was smaller and the quiescent potential less negative when recorded in single β-cell preparations (Rorsman and Trube, 1986; Santos et al., 1991; Kinard et al., 1999; Bertram et al., 2000). It is conceivable that the action potentials might be damped under the patch-clamp recording because the current leak through the gigaseal (∼10 GΩ) between the patch electrode and the cell membrane is comparable to the whole cell membrane current (input resistance of ∼10–30 GΩ). The membrane capacitance of ∼6 pF of small β cells is also in the same order as the floating capacitance of the electrode tip. Moreover, recovery from dissociation injury might be incomplete in culture medium, or action potential generation might be depressed at room temperature or by the rundown of ICaV. It should be noted that mouse or human β cells contain a relatively high density of ICaV ranging over 6 to 11.4 pA pF−1, with an apparent reversal potential of ∼50 mV at physiological [Ca2+]o (see Table I for references). The current densities of ICaV guarantee a fast rising phase and a full size of the action potential in an intact cell before patch-clamp recording. In the present model, the action potential peak was shifted to positive potentials when the membrane K+ conductance was partially blocked, in agreement with experiments (Atwater et al., 1979; Santos and Rojas, 1989; Rorsman et al., 1992; Houamed et al., 2010). In addition to the action potential amplitude, it should also be noted that the burst duration in single-cell preparations might be affected by the leak conductance and floating capacitance during patch recordings. The difference between electrical activities in single cells and those of islets might be caused by the above recording artifacts.
Further considerations and limitations of the study
The VL diagram in Fig. 5 indicated prominently large contributions of ICaV during the whole interburst period at both 8 and 16 mM [G]. These contributions are mainly attributable to the increase or decrease in dCaV, which is a pure voltage-dependent gate of ICaV. From the viewpoint that burst–interburst rhythm is principally generated by slow changes in cytosolic substrate concentrations, the role of dCaV is to magnify changes in Vm in the same direction as those induced by other membrane currents under the influence of cytosolic factors. Namely, at 8 mM [G], the slow depolarization induced by changes in IKATP and INaCa increases dCaV, which results in further depolarization. In the early half of the interburst period at 16 mM [G], dCaV is decreased as a result of the negative shift of Vm, which is primarily induced by a decrease in INaCa or ITRPM via the progressive decay of [Ca2+]i. During the late phase, the gradual positive shift in Vm induced by a decrease in IPMCA or INaK increases dCaV to enhance the depolarization. If these secondary contributions of ICaV are excluded from comparison of the membrane currents, the VL diagram indicates that IKATP and INaCa at 8 mM [G], and INaCa, IPMCA, ITRPM, and INaK at 16 mM [G], play major roles in converting variations in the slow cytosolic factors into the Vm change in our model.
The effects of thapsigargin on [Ca2+]i have been examined in several experiments using islet preparations, because blocking the ER might provide important clues as to the role of Ca2+ buffering by the ER in the bursting rhythm. Unfortunately, there have been no experimental data showing the effects of thapsigargin on the electrical activity or on Ca2+ fluctuations in isolated individual β cells. In our single-cell model, the simulation of applying thapsigarin (Fig. 3) was consistent with the accelerated rhythm of the Ca2+ fluctuations recorded in several experimental observations in pancreatic islets (Miura et al., 1997; Gilon et al., 1999; Fridlyand et al., 2003), provided that the rhythm of Ca2+ transient reflects the electrical bursting activity even in the islet preparations. However, it should be noted that in other islet studies, thapsigargin induced a sustained increase in [Ca2+]i (Worley et al., 1994b; Gilon et al., 1999; Kanno et al., 2002), accompanied by a continuous firing of action potentials (Worley et al., 1994b). Furthermore, both of these responses have been observed in the same experimental study (Miura et al., 1997). At present, it might be speculated that our cell model represents only a given population of β cells, or that the cell-to-cell electrical coupling among different populations of cell types within the islet can produce different patterns in the thapsigargin response.
Although the current system was much improved compared with the previous β-cell models by adding individual current components in the molecular level, further refinement of the formulation will be necessary according to new experimental data in future, especially in respect to temperature effects on channel kinetics. Furthermore, compared with the electrophysiological formulations, the description of energy metabolism is quite simplified in our β-cell model. Therefore, it is beyond the scope of this study to reproduce an ultraslow bursting rhythm with periods >5 min, which have been suspected to be of metabolic origin (Henquin et al., 1982; Bertram et al., 2004). Moreover, the effect of [Ca2+]i on ATP production was not considered in our model. Keizer and Magnus (1989) assumed that Ca2+ entry into the mitochondria depolarized the matrix membrane to inhibit ATP production. However, the opposite effect has also been proposed, namely that an increase in [Ca2+]i might facilitate ATP production by activating dehydrogenases within the TCA cycle (Cortassa et al., 2003). The net effect of [Ca2+]i on ATP production is not quantitatively known at present. Therefore, it will be important to include more precise models for glycolysis (Smolen, 1995; Bertram et al., 2004), TCA cycle, and oxidative phosphorylation (Dzbek and Korzeniewski, 2008) in future formulations.
We thank Professor T. Powell for fruitful discussion and for improving the English of this paper.
This work was supported by the Biomedical Cluster Kansai project; a Grant-in-Aid (22590216 to C.Y. Cha and 22390039 to A. Noma) from the Ministry of Education, Culture, Sports, Science and Technology of Japan; and the Ritsumeikan-Global Innovation Research Organization at Ritsumeikan University.
Lawrence G. Palmer served as editor.
C.Y. Cha and Y. Nakamura contributed equally to this paper.