We have developed a detailed mathematical model of ionic flux in β-cells that includes the most essential channels and pumps in the plasma membrane. This model is coupled to equations describing Ca2+, inositol 1,4,5-trisphosphate (IP3), ATP, and Na+ homeostasis, including the uptake and release of Ca2+ by the endoplasmic reticulum (ER). In our model, metabolically derived ATP activates inward Ca2+ flux by regulation of ATP-sensitive K+ channels and depolarization of the plasma membrane. Results from the simulations support the hypothesis that intracellular Na+ and Ca2+ in the ER can be the main variables driving both fast (2–7 osc/min) and slow intracellular Ca2+ concentration oscillations (0.3–0.9 osc/min) and that the effect of IP3 on Ca2+ leak from the ER contributes to the pattern of slow calcium oscillations. Simulations also show that filling the ER Ca2+ stores leads to faster electrical bursting and Ca2+ oscillations. Specific Ca2+ oscillations in isolated β-cell lines can also be simulated.
- Ca2+ oscillations
- endoplasmic reticulum
- pancreatic islets
- theoretical model
a rise in the intracellular free calcium concentration ([Ca2+]i) is a key signal in the initiation of insulin secretion from the pancreatic β-cell. This increase principally results from calcium influx through plasma membrane (PM) Ca2+ channels, which open in response to secretagogues, primarily glucose. The metabolism of glucose through glycolysis and the tricarboxylic acid cycle leads to an increase in the cytoplasmic ATP-to-ADP (ATP/ADP) ratio. This causes closure of ATP-sensitive K+ (KATP) channels followed by depolarization of the β-cell membrane to the threshold potential where Ca2+ channels open, initiating Ca2+ influx (4). These events underlie glucose-induced electrical activity that, in pancreatic islets, consists of Ca2+-dependent action potentials.
There is extensive literature describing β-cell electrical activity and its relationship to [Ca2+]i in intact islets of Langerhans, isolated islet cells, and insulinoma cell lines. Most of the work has been carried out using mouse islets, with some studies using islets from rat, hamster, human, and other species.
Mouse pancreatic β-cells exhibit complex and cyclic spike-burst activity in response to a rise in extracellular glucose concentration. The bursts consist of a depolarized phase of Ca2+-carrying action potentials alternating with a silent phase of repolarization, resulting in oscillations in intracellular Ca2+, which can drive pulses of insulin secretion (28, 37).
The only stimulus required for a complex cyclic spike-burst activity and corresponding [Ca2+]i oscillations in islets and β-cell clusters is elevation of glucose to levels above 5 and less than ∼20 mM. Intermediate glucose concentrations induce two main types of oscillations in mouse pancreatic islets: fast, where the period ranges from 10 to 30 s, and slow, with periods of several minutes (37, 54, 83). Single mouse β-cells can also respond to glucose stimulation with regular oscillations (37).
We have previously studied slow and fast [Ca2+]i oscillations in islets in response to a variety of conditions (70, 73; unpublished observations). We have also previously reported that a stable, transgenically derived murine insulinoma cell line (βTC3-neo) responds to glucose with slow, large amplitude [Ca2+]i oscillations but only in the presence of 10–20 mM tetraethylammonium (TEA), a blocker of K+ channels (74). We have utilized this cell line to characterize glucose-stimulated oscillatory activity (74).
However, the precise interpretation of previous results is limited due to the numerous channels and pumps in β-cells that work concurrently, and identification of physiologically slow variables that drive oscillations remains unclear. To clarify these complex experimental results, we used a mathematical modeling approach. Our goals, then, are twofold: to develop a model for β-cell ion homeostasis, including the bursts and [Ca2+]i oscillations that can simulate cellular behavior, and to explain on this basis the experimental data for single cells and islets.
Several mathematical approaches in the literature have provided quantitative estimates of glucose-induced insulin secretion with corresponding descriptions of glucose transport, metabolism, and ion regulation (most recent are Refs. 6, 7, 13, 14, 26, 30, 62, 65, and 78). Many phenomena were successfully explained using these models (see DISCUSSION). However, most models described one specific phenomenon and therefore included a very restricted set of channels and pumps, so that it is difficult to apply these models to another aim. In addition, newly described experimental results and signaling molecules should be considered.
For this reason, we developed a new model to consider an extended variety of ionic channels and pumps as well as endoplasmic reticulum (ER) calcium sequestration mechanisms that have been identified in β-cells. We simulated whole cell electrical activity and [Ca2+]i, free calcium in the ER ([Ca2+]ER), intracellular Na+ ([Na+]i), cytosolic ATP ([ATP]i), and inositol 1,4,5-trisphosphate ([IP3]i) concentrations by use of parameters derived mainly for mouse β-cells. However, in this article, we do not make special consideration of metabolic processes or insulin secretion.
To our knowledge, this article represents the first attempt to include in one model the key ion channels and pumps on both the PM and the ER membrane. Several detailed and comprehensive models have been developed for cardiac muscle and other cell types (67, 71). Here, we employ this approach to investigate the possible role of K+ channels, Na+, and IP3 in regulating β-cell Ca2+ oscillations.
MATERIALS AND METHODS
Islet Isolation and [Ca2+]i Measurement
Isolation of rodent islets. Islets of Langerhans were isolated from the pancreata of 8- to 10-wk-old C57BL/6J mice by collagenase digestion and discontinuous Ficoll gradient methods described previously (85).
Measurement of [Ca2+]i. Islets were loaded with fura 2 for 25 min at 37°C in growth medium (RPMI 1640–10% fetal calf serum, penicillin, and streptomycin) supplemented with 5 mM acetoxymethyl ester of fura 2 (Molecular Probes). [Ca2+]i was estimated as described elsewhere (74). Dual-wavelength digitized video fluorescent microscopy with fura 2 in single islets was performed using an intensified charge-coupled device (Hammatsu C2400) and Metafluor imaging software (Universal Imaging).
Our model of the β-cell combines a parallel conductance membrane model and an inside fluid compartmental model (Fig. 1). The compartmental model describes the time rate of changes in [Ca2+]i, [Ca2+]ER, [Na+]i, and [ATP]i, and [IP3]i. The extracellular space is assumed to have a relatively large volume, so that ionic concentrations of Ca2+, Na+, and K+ there are assumed to be constant.
The corresponding equations for particular channels or pumps that we found in the literature were used even though they were in some cases obtained from cells other than pancreatic β-cells. However, we have considered the origin of the equations employed and their correspondence to biological processes. Model coefficients from the literature were adjusted to simulate the corresponding experimental data as indicated. The basic set of evaluated coefficients is shown in Table 1, denoted as “standard.”
Plasma membrane currents. The differential equation describing time-dependent changes in the plasma membrane potential (V) is the current balance equation (41) 1 where Cm is the total membrane capacitance. Model currents include a voltage-dependent Ca2+ current (IVCa), a PM calcium pump current (ICa,pump), Na+/Ca2+ exchange current (INa,Ca), Ca2+ release-activated nonselective cation current (ICRAN), inward Na+ currents (INa), a Na+/K+ pump current (INa,K), a delayed rectifier K+ current (IKDr), voltage-independent small conductance Ca2+-activated K+ current (IKCa), and the KATP current (IKATP) (Fig. 1).
Ionic flux across the membrane, Jxi, is related to the ionic current by the expression 2 where F and zx, are Faraday's constant and ionic charge, respectively.
General current equation. The ionic current across channels may be given by the empirical equation (41, p.18) 3 where 4 5 and where gXi is the whole cell conductance, VX is the equilibrium potential, R is the universal gas constant, T is absolute temperature, and Xo and Xi are the extra- and intracellular concentrations, respectively.
Whole cell conductance may depend on ion concentrations and voltage. These dependencies will be specified for each conductance as they are discussed.
IVCa. The whole cell Ca2+ current in mouse β-cells flows principally through voltage-activated Ca2+ channels (29, 80). This current increases when the membrane is depolarized. Equation 3 was used for this current. The current/voltage relationship of this current [pVCa(V)] can be described by a Boltzmann-type activation curve (7, 9). Then 6 where 7 and where gmVCa is the maximum whole cell conductance, VCah is the half-activation potential, and KCah is the slope at half-maximal potential.
At physiological Ca2+ concentrations, the measured VCah varies from -3.8 mV (80) and -16.7 mV (45) to -19 mV (9), and KCah extends from 6.6 mV (45) and 8.4 mV (80) to 13 mV (9). In our model, these coefficients were fitted to be inside these regions (Table 1).
According to Göpel et al. (33) the mean integrated Ca2+ current observed during a 100-ms depolarization to -20 mV in the β-cell in situ was 7.7 pC at [Ca2+]o = 2.6 mM. This result can be obtained from Eq. 6 of our model, with gmVCa = 995 pS and [Ca2+]i = 0.05 μM. We set gmVCa = 770 pS as the standard condition for our model (Table 1). This value is close to that estimated above. The current/voltage relationship calculated using Eq. 6 is in good correlation with data for β-cells in situ (33).
PM Ca2+ pumps (ICa,pump). PM Ca2+ pumps are now well recognized as a primary system for the specific expulsion of Ca2+ from β-cells (87). However, their kinetic properties are not well studied in these cells. For this reason, the general properties of PM Ca2+ pumps were used for the calculations.
It is now generally established that the Ca2+/ATP stoichiometry of the PM Ca2+ pump is 1 (10, 66) and that the reconstituted PM pumps are capable of establishing a membrane potential while operating with H+/Ca2+ = 1 (40). PM and ER Ca2+ pumps do not have the sensitivity to ATP at the millimolar ATP concentrations that exist in physiological conditions (10). The Km is <0.5 μM Ca2+ for activated PM pumps (10, 66). On the basis of these data, the Ca2+ current through PM Ca2+ pumps was expressed in our model as (see 62) 8 where PmCap is a maximum current and KCap is the value for the half-activation calcium concentration. KCap was set to 0.1 μM, which is within the range of experimental observations for other preparations. PmCap (Table 1) is correspondingly above the evaluated lower limit (830 fA) that we obtained (see RESULTS).
Na+/Ca2+ exchanger (INa,Ca). The lack of specific inhibitors hinders the experimental assessment of the role of the Na+/Ca2+ exchanger. However, a detailed investigation of the β-cell electrogenic 3 Na+/Ca2+ exchanger was recently reported (25), and the exchanger current (INa,Ca) was taken as follows from this work 9 where 10 and where gNaCa is the maximum whole cell conductance, and KNaCa is affinity constant.
One can calculate from Gall et al. (25) that gNaCa = 271 pS (as Cm·ğNa/Ca from their Eq. 2), and we used this value for gNaCa in our model.
Similarly, those authors found KNaCa to be 1.5 μM (25). However, an insignificant INa,Ca and a correspondingly very small influence of the Na+/Ca2+ exchanger on the model solutions were obtained at KNaCa = 1.5 μM. This result does not correspond to the data involving the effect of this exchanger on glucose-induced electrical activity in intact pancreatic islets (25). This apparent contradiction is probably due to the use of low ATP concentration in the perforated patch whole cell experiments in Ref. 25; therefore, KNaCa (1.5 μM) should only be considered as an upper limit. Data from different cell types show that ATP induces a dramatic increase in the intracellular Ca2+ affinity for the Na+/Ca2+ exchanger (20). Correspondingly, a lower value for KNaCa (0.75 μM) was employed in our model to fit glucose-induced β-cell oscillation patterns at an increased ATP content.
ICRAN. ICRAN is a nonselective cation current whose conductance is regulated by the ER Ca2+ content. In mouse β-cells, it could be activated either indirectly by ER Ca2+ store depletion or directly by maitotoxin (84, 88, 89). However, the mechanism for coupling ER Ca2+ store depletion with CRAN remains an unresolved question.
In our model, only Na+ can penetrate into the cell via this channel, reflecting experimental results under physiological concentrations of cations (52, 84, 88). The current/voltage relationship of this current [pCRAN(V)] was roughly linear, with a reversal potential close to 0 mV (52, 84, 88). In previous models, it was suggested that the regulation of ICRAN depends on [Ca2+]ER (8, 14, 62). We take this dependence following Ref. 62. Equation 3 can be written for this Na+ inward current 11 where 12 13 and gmCRAN is the maximum whole cell conductance, fCRAN is the voltage-independent part of the current, KCar is the half-activation [Ca2+]ER level, and VCRAN is the reversal potential.
INa. The general equation for INa is derived from Eq. 3. The relationship between conditioning voltage and relative current amplitude for tetrodotoxin-blockable component of the INa [pNa(V)] was characterized (33) as the Boltzmann equation. Then 14 where 15 The maximum value of whole cell conductance (gmNa) was evaluated from the data (33), where it was found that INa = -392 pA at V = -140 mV. This INa value can be obtained from Eq. 14, with gmNa as 1,200 pS when [Na+]i is close to 0 and [Na+]o = 140 mM. (We used [Na+]i = 0.1 mM for this calculation). This value was used for standard conditions (Table 1).
Na+-K+ active transport (INa,K). Electrogenic Na+-K+-ATPase extrudes three Na+ ions in exchange for two K+ ions for each molecule of ATP hydrolyzed, generating a net outward flow of cations through the β-cell PM (68). However, its kinetic properties have not been studied in β-cells. For this reason, we employed the general model (12) as was done previously for β-cells (65) (see APPENDIX A).
Equation 3 was used; however, the formulation (26, 78) is employed for whole cell conductance. 16 where g′mKDr is the maximal whole cell conductance, and n is the general gating variable, which is determined in APPENDIX B.
Family of KCa channels (IKCa). The voltage-independent small-conductance Ca2+-activated K+ channels (SK channels) are also expressed in β-cells (unpublished observations). It is quite plausible that these channels play an important role in [Ca2+]i oscillations (32).
Equation 3 was used for this channel. Ca2+ activates IKCa currents [f(Ca2+)] in a sigmoidal fashion (49, 53). However, there are no data on the exact dependence in β-cells. The Hill equation was used to model this relationship, where the Hill coefficient equals 3 (14) or 5 (62). We also use a Hill equation, but with a Hill coefficient equal to 4 17 where 18 and KKCa is affinity constant.
KKCa varies in the 0.05–0.9 μM range (Ref. 41, p. 144). However, SK3 channels may be the predominant SK family member in β-cells (unpublished observations), and because KKCa = 0.1 μM was found for SK3 channels by Carignani et al. (11), this value was used in our model.
It is unlikely that the large-conductance Ca2+-activated K+ channel (BK channel) encoded by the slo gene plays an important role in the generation of oscillatory activity in β-cells (50), although several previous models of β-cell bursting have included this channel (78). We therefore did not include it in our considerations.
KATP channels (IKATP). Equation 3 was used for this channel. However, we adopt here a kinetic model for the value of whole cell conductance (see Refs. 43, 58, 65). Then 19 where g′mKATP is the maximal value of the whole cell conductance, and OKATP is the fraction of channels open, which is considered in APPENDIX C.
The value g′mKATP was evaluated to be 24 nS for the standard coefficient, which is comparable to the value used in Refs. 59 and 65. Then it can be calculated from Eqs. 19 and C1 (APPENDIX C) that the value of the whole cell conductance (g′mKATP OKATP) is 0.94 nS at [ATP]i = 0. This value falls in the range between 0.2 and 3.0 nS that was measured for the maximum KATP conductance reported for single mouse β-cells (48).
Fluid Compartmental Model
The material balance equation for calcium depends on three processes: entry, extrusion, and buffering (Fig. 1).
Calcium enters the β-cells primarily through voltage-activated Ca2+ channels by diffusion along an inwardly directed electrochemical gradient (IVCa). The maintenance of this ionic gradient depends on a Ca2+-extruding mechanism in the PM and Ca2+ sequestration by intracellular organelles. At the PM, three processes are involved in transporting Ca2+ out of the cell: a Ca2+ pump, a Na+/Ca2+ exchanger, and removal of Ca2+ sequestrated in insulin granules by exocytosis. In addition, both the ER and mitochondria can accumulate Ca2+ via pumps. Even though Ca2+ is critical for mitochondrial function, for the present we do not include mitochondrial Ca2+ stores, since it appears that both the volume of mitochondria and its Ca2+ concentration are small relative to the ER (2, 18).
Sarco(endo)plasmic reticulum Ca2+-ATPase pump. Uptake of Ca2+ into the ER is mediated in β-cells by sarco(endo)plasmic reticulum Ca2+-ATPases (SERCAs) with an unusually low value for the half-activation calcium concentration (≤0.4 μM) (27). The Ca2+/ATP stoichiometry of the SERCA is 2 (10, 66). The equation for SERCA was written in the same form as the Ca2+ PM pump (Eq. 8; see also Ref. 62) 20 where Jer,p is the flux into the ER through SERCA pumps per cytosol volume, KCarp is the half-maximal pump activity, and PCaER is the maximum rate of pumping. The lower limit for this parameter was evaluated (see RESULTS) as PCaER = 0.026 μM/ms for the data from Ref. 25. The evaluated parameter for standard conditions (Table 1) is correspondingly above this lower limit. KCarp was taken as 0.5 μM, a value close to experimental data (27).
IP3 metabolism. In the pancreatic β-cell, there is clear evidence that IP3 is an important cellular messenger inducing Ca2+ mobilization from intracellular stores by binding to specific receptors located on intracellular Ca2+ stores and is the main calcium release channel from the ER in β-cell tissues (38). However, we could not find any mathematical models of β-cell regulation that used IP3 as an independent variable.
Different mechanisms seem to be available in the β-cells for stimulating IP3 production; however, only the modulation of inositol lipid-specific phospholipase activity by changes in Ca2+ is well studied (3). For this reason, the [Ca2+]i dependence was only modeled as the Hill function with a Hill coefficient of 2 (3).
The subsequent kinetics of IP3 in β-cells is unknown. Therefore, we assume simply that IP3 is converted to inactive metabolites at a rate proportional to its concentration. Then synthesis and degradation of IP3 are described by the equation 21 where kIP is the rate constant of IP3 production, kdIP is the rate constant of IP3 degradation, and KIPCa is the half-activation of the IP3 production by calcium concentration.
The time constants kIP and kdIP (Table 1) were picked to correctly reproduce the experimentally observed time course of the average [Ca2+]i during the slow [Ca2+]i oscillations (see RESULTS). KIPCa is taken as 0.4 μM (36). The fitted parameter kdIP (0.04 s-1, Table 1; half-life is 17.3 s) is inside the region of a half-life from 1 to 30 s that has been found for different cells with a radius of ∼5 μm and above (79), similar to that of β-cells.
CaER2+ mobilization. A Ca2+ leak current (Jout) from the ER per whole cell was taken according to Mears et al. (62) 22 where Pleak is the Ca2+ leak permeability out of the ER, PIP3 is the maximum permeability through open IP3-activated channels, and O∞ is the fraction of open IP3-activated Ca2+ channels.
The type III IP3 receptor (IP3R) is the main calcium release channel in the β-cell ER (38). Apparently, in normal physiological conditions, the O∞ does not depend on the change of ATP concentration (57). The initial dependence of the IP3R maximum open probability vs. cytosolic Ca2+ could be well fitted to a Michaelis-Menten-type saturating equation with an affinity of 0.077 μM for [Ca2+]i (60).
However, there are conflicting data concerning IP3 action. According to Hagar and Ehrlich (39), the Hill equation with a Hill coefficient close to 2 was found for the effect of IP3 on O∞ in vitro. On the other hand, the Hill number for low degrees of saturation of the IP3R was estimated at ∼3.5 in vivo for rat basophilic leukemia cells (63). On the basis of these data, the Hill coefficient is set equal to 3. Then 23 where KRCa is the affinity constant of the [Ca2+]i to IP3R, and KIP3 is the constant for IP3 binding to a channel. KIP3 is 3.2 μM (39).
Ca2+ and Na+ dynamics. The cytosolic and ER calcium concentrations and Na+ cytosolic concentration are determined by the ion fluxes across the plasma and ER membranes. However, calcium concentrations are strongly buffered in cells (77). This can be modeled using special coefficients for the fraction of free Ca2+ in the cytoplasm and ER (14, 62, 71). The Ca2+ leak from cells by β-granule exocytosis was modeled previously (14) as a rate proportional to [Ca2+]i. We employed a similar approach.
On the basis of the foregoing consideration, the equations for Ca2+ and Na+ concentrations can be written as 24 25 26 where fi and fer are the fractions of free Ca2+ in cytoplasm and ER, Ver and Vi are the effective volumes of the ER and cytosolic compartments, and ksg is a coefficient of the sequestration rate of Ca2+i by the secretory granules.
It has been proposed that the Ca2+-binding capacity inside the ER is lower than in the cytoplasm (69). This can increase the fraction of free Ca2+ in ER; therefore, we used a value of fer threefold higher than fi (Table 2).
ATP homeostasis. ATP concentration in the cytosol increases somewhat with an increase in glucose supply (2). However, the rate of ATP production and utilization is a complex function of the concentrations of ATPi, ADPi, Ca2+, glucose, and numerous other factors. For simplicity, we have chosen a first-order reaction to express the rate of ATPi production from ADPi.
Clearly, there is some basal level of ATPi consumption in a cell at low glucose and [Ca2+]i. There is also ATP consumption by the PM and ER Ca2+ pumps and by the electrogenic Na+-K+-ATPase in our model. However, there is considerable evidence that other ATP consumption processes are accelerated by an increase in [Ca2+]i in β-cells (2, 23, 76). For this reason, an additional term for ATP consumption was introduced, which depends on an increase in [Ca2+]i. Then 27 where 28 kADP is the rate constant of ATPi production, kATP,Ca is the rate constant of ATPi consumption that accelerates as [Ca2+]i increases, and kATP is the permanent rate constant of ATPi consumption.
The general concentration of intracellular nucleotides ([A]o) was assumed to be constant. The measured concentration of adenine nucleotides in pancreatic islets varies within the limits of 2–10 mM (23). We have taken [A]o to be inside this range (4 mM).
It is important to find a range of parameters for ATP dynamics. According to the calculations in Ref. 23, the rate of ATPi production in β-cells may run to ∼0.5 μmol ATP·s-1·g dry wt-1 at low glucose concentrations and 1 μmol ATP·s-1·g dry wt-1 at high glucose concentrations. Matschinsky et al. (61) found the corresponding level of adenine nucleotides to be ∼20 μmol/g dry wt. Then the mean rate of ATPi production per second can be estimated roughly as a part of the general quantity of intracellular nucleotides per unit of volume, that is, as 0.025 or 0.05 [Ao]/s (or 0.1 and 0.2 mM/s) at low and high glucose concentrations, respectively.
The constants kATP,Ca and kATP (Eq. 27 and Table 1) were fitted to simulate both the observed pattern of [Ca2+]i oscillations and the evaluated rate of ATPi production in β-cells above. The simulation of low glucose concentration (kADP = 0.03 s-1) leads in our model to a steady-state [ADP]i of 3.068 mM (Table 3) and to a corresponding rate of ATPi production (kADP [ADP]i) of 0.092 mM/s. The simulation of slow [Ca2+]i oscillations at high glucose levels (kADP = 0.37 s-1; Fig. 3) yields a mean [ADP]i of 0.7 mM and to the rate of ATPi production as 0.26 mM/s. This correlates well with the values evaluated above.
Computational aspects. The spread of current between electrically coupled β-cells most probably contributes to the synchronization of the electrical activity and [Ca2+]i oscillations throughout an islet (33, 37). For this reason, we consider islets as an assemblage of cells with similar properties and perform the computer simulation only for some mean individual cell.
The complete system consists of seven state variables: differential equations describing the time rate of change in PM potential (V), Ca2+, Na+, ATP, and IP3 concentrations in cytoplasm, calcium concentration in ER, and the differential equations characterizing the voltage-dependent gating variable (n) for the delayed rectifier K+ channels (Eqs. 1, 21, 24–27, B1). The units used in the model are time in milliseconds (ms), voltage in millivolts (mV), current in femtoamperes (fA), concentration in micromoles per liter (μM), conductance in picosiemens (pS), capacitance in femtofarads (fF), and temperature in degrees Kelvin (°K).
The model and cell parameters are given in Tables 1, 2, 3. These tables and the appendixes contain all the information necessary to carry out the simulations presented in this paper. However, the rate constant for ATPi production (kADP) is represented in every figure legend.
Simulations were performed by forward integration of the coupled system of differential equations. We used two integrated software environments: “Content” for the IBM-compatible personal computer using an implicit fourth-order Runge-Kutta or Rodau IIA method (15) and “Virtual Cell” accessible via the Internet (55).
This model is available for direct simulation on the website “Virtual Cell” (www.nrcam.uchc.edu) in “MathModel Database” on the “math workspace” in the library “Fridlyand” with the name “Chicago.1”.
Validation of the Model
It was possible to evaluate the values of the rates of Ca2+ pumps from the simulation of the data from Ref. 25, where the [Ca2+]i transients were induced by a 2-s depolarization. Gall et al. (25) found that the time constant for the decrease in [Ca2+]i is ∼1.8 s when Ca2+-ATP pumps on the PM and ER were activated and 4.6 s when Ca2+-ATP pumps on the ER were inhibited by thapsigargin (Tg). This was determined under conditions where the Na+/Ca2+ exchanger was inactive, i.e., gNaCa = 0, granule exocytosis was lacking (ksg = 0), and IP3R was closed (PIP3 = 0). Using Eqs. 9 and 21, we found that these data can be well simulated at PCaER = 0.026 μM/ms and PmCap = 830 fA. However, these values should be considered as a lower limit, whereas the experiments in Ref. 25 were made in conditions of low ATP concentration, and it is probable that Ca2+ pumps were in an inactivated form. Correspondingly, higher values were fitted in our model (see Table 1).
A number of parameters were taken from the literature, such as channel and pump kinetics (see MATERIALS AND METHODS). However, because these values were sometimes derived from experiments with different cell lines and specific experimental conditions, these parameters were constrained to be within a range consistent with our own system. Other parameters were fitted to simulate the pattern of fast and slow bursting of electrical activity and [Ca2+]i oscillations in β-cells. The parameters that yield these oscillations are the standard values given in Table 1. A computer simulation at a low ATP production rate and the coefficients from Tables 1 and 2 were used to obtain the steady-state data corresponding to low glucose concentrations (Table 3).
Slow Ca2+ oscillations. Slow Ca2+ oscillations occur spontaneously during glucose stimulation (Fig. 2). They are characterized by a descending phase in two components. A rapid decrease in [Ca2+]i is followed by a slower phase. Similar data were obtained by many investigators for slow oscillations (27, 37, 64). It was also shown that a rapid decrease in [Ca2+]i coincided with the closure of voltage-dependent Ca2+ channels (27).
A typical computer simulation of slow oscillations with the patterns of membrane potential and corresponding concentrations is illustrated in Fig. 3, left. It was obtained by introducing a step increase in the rate constant of ATPi production from the resting conditions (Table 3) to an intermediate value. A detailed pattern of the concentration changes and corresponding currents is shown in Fig. 5A for a single oscillation.
As can be seen in Fig. 3, the acceleration of ATPi production increases the [ATP]i/[ADP]i ratio after a step increase of kADP. This closes KATP channels in the PM. This closure decreases the IKATP that depolarizes the PM, leading to the opening of voltage-dependent Ca2+ channels, the stimulation of Ca2+ influx, and, eventually, a rise in [Ca2+]i. Increasing the IVCa makes a contribution to depolarization of the PM by itself. The beginning of electrical spiking is generated at threshold PM potential to form bursts and leads to [Ca2+]i increase (Fig. 3B, arrow 1). Continuous spiking follows up to a point at which the repolarizing processes lead to decreased PM potential. This in turn leads to the closure of voltage-gated Ca2+ channels and the continuous spiking terminates (Fig. 3B, arrow 2). In the silent phase, slow depolarization of the PM begins until all processes repeat, leading to continuous oscillations. Dissection of the burst mechanism and the effects of Tg will be discussed in Model Sensitivity to Slow Variables.
Biphasic responses and fast Ca2+ oscillations. Fast oscillations usually show a more complex picture. In this case, the electrical response of pancreatic islets to step increases in glucose concentration is often biphasic, consisting of a prolonged depolarization with action potentials (Phase 1) followed by fast bursts of membrane potential and Ca2+ oscillations. The initial tonic electrical activity is accompanied by a sustained elevation of [Ca2+]i (8, 72, 89).
It was suggested by model simulations and experimental results that the Phase 1 response may result from the combined influences of depolarizing KATP channel closure and ICRAN deactivation. The Ca2+ fills the ER during Phase 1, deactivating ICRAN and repolarizing the PM, allowing steady-state bursting to commence (8, 62). We are also able to simulate these results on the basis of our model when Ca2+ fills the ER to considerably elevated concentrations. This increase of [Ca2+]ER can occur if Ca2+ output from the ER is restricted, for example, by closing IP3R channels. This can be achieved in our model, for example, by decreasing the [IP3]i as a consequence of decreased IP3 production. The typical Phase 1 response followed by spike-burst behavior of the membrane potential and the fast oscillation patterns of ions and [ATP]i are illustrated in Fig. 4. It was simulated at a decreased IP3 synthesis rate constant and an intermediate value of the rate constant of ATPi production.
The same results can be simply obtained by a partial closing of the IP3R channels in place of the decreased IP3 production (not shown). Our model also shows that the duration of Phase 1 decreases with increased initial [Ca2+]ER at t = 0 (not shown; see Ref. 62).
Simulating increased ATPi production leads to increased oscillation frequency, progressing to continuous spiking with a further increase in the ATPi production rate constant (Fig. 4). This is in reasonably good agreement with the experimental data (13, 37).
Note that fast [Ca2+]i oscillations simulated by our model (Fig. 4) have an appearance very similar to those observed in pancreatic β-cells (8, 72). They also resemble the fast oscillations that were obtained in preceding models (8, 14, 62). This is not surprising, because our simulation of fast oscillations is close to those employed in these models.
Model Sensitivity to Slow Variables
Four parameters can change slowly in our model ([ATP]i, [Ca2+]ER, [IP3]i, and [Na+]i), and their variations can lead to a slow membrane repolarization during a depolarized phase and to a depolarization in a silent phase. To determine the role of different processes and parameters, we performed an analysis of model sensitivity to slow variables to fix some parameters at their mean level for some kind of oscillations. It was found that the slow oscillations take place when only [Na+]i ([ATP]i, [Ca2+]ER, and [IP3]i were constant) can change (Fig. 5D), and the characteristic fast oscillations occur when only [Ca2+]ER ([ATP]i, [Na+]i, and [IP3]i were constant) can vary (Fig. 5C). Consideration of these two cases simplifies the general analysis.
The behavior of [Ca2+]i, V, [Ca2+]ER, [IP3]i, and [Na+]i and the most essential current components is shown in Fig. 5, A and B, for single oscillations of slow and fast types. Because of the large fluctuations of ICRAN, the voltage-independent part of this current (fCRAN, Eq. 12) is also represented to facilitate the analysis (broken lines in Fig. 5.9). Corresponding concentrations and currents are also represented for cases when only [Ca2+]ER (Fig. 5C) or [Na+]i (Fig. 5D) can change. The contribution of other currents was either insignificant (INa) or does not change considerably between cases (ICa,pump, IKDr, IKCa, IKATP).
ATP sensitivity. Activation of the Na+/K+ pump with increased [Na+]i concentration during the depolarized phase accelerates ATP consumption in our model. Total Ca2+ pump rates and Ca2+-dependent ATPi consumption also increase during the depolarized phase with [Ca2+]i rising (see Eq. 27). The acceleration of ATPi consumption leads to a decreased [ATP]i/[ADP]i ratio (Figs. 3 and 4). It is therefore possible that cyclic changes in the [ATP]i/[ADP]i ratio control the cyclic closure and reopening of KATP channels. This has been proposed as a mechanism underlying the oscillatory behavior of β-cells (1, 65, 82).
We evaluated this possibility and found that fixing [ATP]i at a mean level yields no noticeable effect on the oscillatory pattern shown in Fig. 5, A and B. This can be explained by insignificant changes of the fraction of KATP channels open (APPENDIX C) during the ATPi variations around its mean value (not shown). This suggests that [ATP]i is not the variable that drives slow oscillations in our model. On the other hand, considerable changes in mean [ATP]i and correspondingly in IKATP can lead to a dramatic change in oscillations (see Figs. 3 and 4).
Single-oscillation analyses. It can be seen that spike activity is generated by the alternate activation of delayed rectifying K+ channels and voltage-gated Ca2+ channels in every case. It is analogous to other models and has been previously analyzed (14, 78). The main difference between the models includes the nature of processes that lead to slow depolarization of PM potential during the depolarized phase and to its repolarizing during a silent phase.
When [Ca2+]ER is the only slow parameter ([ATP]i, [Na+]i, and [IP3]i were constants; Fig. 5C), the following mechanism is responsible for oscillations: rising Cai2+ during the depolarized phase leads to increased Ca2+ pumping from the cytoplasm to the ER and to CaER2+ accumulation, with corresponding closure of CRAN channels and decreased ICRAN, that can be observed more clearly from the change of the voltage-independent part of ICRAN (fCRAN, broken line in Fig. 5C.9). This, in turn, leads to PM repolarization. These processes have opposite directions in the silent phase.
Figure 5D shows that, when [Na+]i is the only slow variable ([ATP]i, [Ca2+]ER, and [IP3]i were constant), INa,K is the main current that manages the slow PM potential changes leading to oscillations. The mechanism involves increased [Ca2+]i during a depolarized period, activating the outward Ca2+ flux through Na+/Ca2+ exchangers and a corresponding inward Na+ flux. The resultant increase in [Na+]i leads to the slow increase of outward current through electrogenic Na+/K+ pumps (INa,K; Fig. 5D.7) with corresponding PM repolarization. Falling [Na+]i during a silent phase leads to decreased INa,K and then in turn to PM depolarization.
In the simulation of slow oscillations shown in Fig. 5A, [Na+]i and the current INa,K change as in Fig. 5D, with no change of the voltage-independent part of the current ICRAN (fCRAN). Correspondingly, the mechanism of slow oscillations is the mechanism described above, with [Na+]i as the only slow variable.
In the simulation of fast oscillations shown in Fig. 5B, [Ca2+]ER and fCRAN change as in Fig. 5C, with little change of [Na+]i and the corresponding current INa,K. In this case, the currents INa,K and ICRAN act together. However, ICRAN plays the major part, as the concentrations and currents are similar in Fig. 5, B and C. This is in agreement with the conclusion that [Ca2+]ER drives the fast oscillations but does not need to oscillate to have slow oscillations, whereas [Na+]i gives slow oscillations but does not need to oscillate to have fast oscillations.
However, in the case of fast oscillations, CRAN channels are working under conditions when the ER is filled with Ca2+, and the [Ca2+]ER change occurs close to the half-activation level of [Ca2+]ER for CRAN channels (200 μM; see Eq. 12 and Table 1). This leads to an increased rate of both PM potential depolarization and repolarization, when [Ca2+]ER correspondingly increases and decreases following [Ca2+]i changes.
Shape of slow [Ca2+]i oscillations and role of IP3. It can be seen in our simulations that, at the beginning of the silent period of slow oscillations, CaER2+ is gradually released into the cytosol (Figs. 3 and 5A.3) and the shape of the [Ca2+]i curve is close to what was found by us (Fig. 2) and others (27, 83) in experiments. Such a slow phase is absent at the simulation of fast oscillations (Figs. 3 and 5B.1).
When [IP3]i was fixed at its mean level at the conditions that were used to simulate slow oscillations, there was still the ability to produce some oscillations. However, the slow [Ca2+]i decrease disappears at the beginning of the silent period. This resulting single [Ca2+]i oscillation is shown as a broken line in Fig. 5A.1.
The relationship between [Ca2+]i and [Ca2+]ER needs additional clarification. In our model, a dynamic equilibrium is established between [Ca2+]i and [Ca2+]ER for several seconds when there is no spiking activity; that is, [Ca2+]i and [Ca2+]ER are closely connected. The decrease in [IP3]i leads to a corresponding, slow IP3R channel closure that decreases the rate of Ca2+ flux from the ER. A slow [Ca2+]ER drop leads to a slow decrease in [Ca2+]i at their dynamic equilibrium. When the [IP3]i is fixed, this mechanism does not operate, and the slow [Ca2+]i drop disappears.
No remarkable changes in the fast oscillations pattern (shown in Fig. 5B) were observed if [IP3]i was fixed at its mean level with otherwise identical conditions (not shown). The lack of an effect can be understood in light of the small changes in [IP3]i that occur during fast oscillations (see Figs. 4 and 5B.4).
This suggests that changes in [IP3]i are not the main driving force for the fast and slow oscillations in our model. However, the change of mean [IP3]i can have a considerable effect, including a transformation of slow oscillations to fast at decreasing [IP3]i (compare Figs. 3 and 4). Furthermore, [IP3]i changes can have a considerable effect on the pattern of slow oscillations.
Role of the ER. In islets pretreated with Tg to block the SERCA pump, the frequency was increased, the amplitude of slow [Ca2+]i oscillations was larger than in control islets, and the descending phase of each [Ca2+]i oscillation was fast with no slow descending second phase (Fig. 2; see also Refs. 27, 64). Our model permits a simulation of these experiments (Fig. 3). In this case, the corresponding amplitude of slow [Ca2+]i oscillations was larger, the period of oscillation was shorter, and the descending phase of each [Ca2+]i oscillation was fast with no slow second phase (Fig. 3), all of which are in reasonably good agreement with the experimental data.
In this case, Ca2+ in the ER does not drive the oscillations. However, the remaining [Na+]i-dependent mechanisms are able to create the relevant oscillations, and corresponding fluxes resemble the ones that were observed if the oscillations were driven only by Nai variations, as in Fig. 5D.
A different pattern was found when SERCA pumps were blocked in fast oscillation experiments. Worley et al. (89) observed that, when Tg was added to the medium, the bursting was transformed to continuous spiking, and fast [Ca2+]i oscillations were converted to a tonically elevated Ca2+ level. Bertram et al. (8) found that an addition of Tg was accompanied by a gradual increase in the plateau fraction and in electrical burst frequency.
However, in our general model, blocking of SERCA leads to the pattern of slow oscillations, as shown in Fig. 3, right. For this reason, other possibilities were investigated. When oscillations were simulated only by changes in [Ca2+]ER (by fixing [Na+]i, [IP3]i, and [ATP]i to constant values, as in Fig. 5C), the block of the SERCA pumps increased the frequency of fast [Ca2+]i oscillations, because ICRAN is activated in consequence of [Ca2+]ER decrease. A further increase of SERCA inhibition can then lead to a permanent PM depolarization and to a continuous rise in [Ca2+]i (Fig. 6). This transformation is consistent with the aforementioned experiments for fast oscillations. It is also in agreement with the results of the modeling of Chay (14) and Bertram et al. (8), who had not included [Na+]i changes in their models and had taken into account only Ca2+ and K+ currents and [Ca2+]ER influence on ICRAN.
In other words, the model suggests that, under conditions when fast oscillations occur, something is blocking the [Na+]i-dependent oscillation mechanisms, even though Tg action resulted in a depletion of the CaER store. One possible mechanism is a partial inhibition of the Na+/K+ pump. For example, we found that only [Ca2+]ER-dependent fast oscillations were possible in our model when the maximum activity of the Na+-K+-ATPase was diminished to 25% of basal activity, even though all other coefficients and initial concentrations were as in Fig. 4. In this case, a simulation of Tg action at fast oscillations resembles Fig. 6. However, the physiological differences between fast and slow oscillations need further examination.
Role of [Na+]i. Slow glucose-induced [Na+]i oscillations were found in individual β-cells under conditions known to induce oscillations of [Ca2+]i (34, 35). Partial suppression of the Na+/K+ pump by ouabain resulted in an increased amplitude of slow Na+ and Ca2+ oscillations and decreased their frequency. Slow [Na+]i and [Ca2+]i oscillations were converted to tonically elevated [Na+]i and [Ca2+]i levels following increased ouabain concentration (34, 35).
Our simulation leads to analogous results (Fig. 7). This result can be explained on the basis of our model: the partial suppression of the Na+/K+ pump leads to decreased INa,K, which is responsible for slow PM potential variations. For this reason, to obtain the same value of INa,K variation that leads to a corresponding PM potential change, it was necessary to have larger [Na+]i changes. This extends the time necessary for corresponding [Na+]i accumulation or dissipation. The increase in Na+/K+ pump inhibition leads to collapse of this [Na+]i-connected mechanism and to the disappearance of slow oscillations.
To check the role of [Na+]i in driving slow oscillations, we also froze [Na+]i at its mean level characterized for slow oscillations (7 mM) and made a simulation at the initial conditions by using the coefficients characterized for slow oscillations (as in Fig. 3). No oscillations were found (not shown).
Our model simulations demonstrate the necessity for [Na+]i oscillations to drive the slow Ca2+ oscillations and suggest that [Na+]i is a candidate pacemaker for slow oscillations.
We have previously studied the properties of a clonal β-cell line termed βTC3-neo (74). In contrast to normal mouse islets, this cell line does not usually respond to a step increase in glucose concentration with regular oscillatory increases in [Ca2+]i. Instead, a slow rise in [Ca2+]i occurs with occasional intermittent spikes. However, the exposure of βTC3-neo cells to 20 mM TEA, a blocker of K+ channels (24), permitted the generation by glucose of large, regular, slow oscillations in [Ca2+]i. In the absence of glucose, TEA was without effect (74). However, the molecular differences between βTC3-neo cells and primary β-cells have not been completely studied.
The magnitude of the Ca2+ inward current in β-cells within intact pancreatic islets is considerably more than in individual β-cells maintained in tissue culture [up to twice (33)]. This difference can be expected to have dramatic effects on the generation of β-cell electrical activity and leads to failure of dispersed cells to generate the characteristic oscillatory electrical activity (33, 81). It is reasonable to assume that βTC3-neo cells also have a decreased Ca2+ inward current that can be simulated by a reduced conductance of Ca2+ channels in our model.
Under these conditions, i.e., reduced conductance of Ca2+ channels, simulation of glucose addition does not lead to oscillations. However, a partial block of delayed-rectifier K+ channels, small-conductance Ca2+-activated K+ channels, and KATP channels, which simulate TEA action, leads to large, regular oscillations in our model (Fig. 8). This occurs because the decreased inward Ca2+ current leads to PM repolarization. To restore the possibility of oscillations, a corresponding decrease of K+ outward currents, as by TEA action, is required to allow PM depolarization.
The exposure of βTC3-neo cells to 1 μM nitrendipine, which blocks L-type Ca2+ channels, can reversibly suppress oscillatory activity (Fig. 1C in Ref. 74). In our model, a similar decrease in Ca2+ channel activity also leads to a suppression of oscillatory activity (not shown).
In the absence of glucose, exposure of βTC3-neo cells to 100 μM tolbutamide to block IKATP induced [Ca2+]i oscillations only in the presence of TEA (Fig. 2A in Ref. 74). The same results were obtained with our model (Fig. 9) when an additional block of KATP channels was used to simulate tolbutamide effects. In this case, the simulation of tolbutamide addition on IKATP was equivalent to the simulation of glucose additions in Fig. 8.
Similar experimental results were obtained by Eberhardson and Grapengiesser (21) for single β-cells isolated from mouse islets, in which the [Ca2+]i oscillations could be obtained at low glucose concentration only when both TEA and tolbutamide were added.
We have attempted to integrate plasma membrane and intracellular events in modeling the behavior of pancreatic β-cells by building on several previous models. The model includes consideration of ionic flux, metabolism, and regulation by intracellular messengers. The complexity of the system necessarily left us with some flexibility in assigning values of the parameters. Therefore, it would be unreasonable to expect the model to accurately reproduce all previously observed experimental results. However, it is necessary to verify that the model is able to reproduce the general features of β-cell or islet oscillation patterns.
Islet Oscillatory Behavior
The majority of previous models emphasized the features of fast-type oscillations, and detailed theoretical considerations of this kind of oscillations have been made (14, 26, 62, 65, 78). We were also able to simulate these fast oscillations (Figs. 4 and 5B). Slow oscillations most likely constitute a physiological oscillatory pattern in β-cells, and they may constitute the framework for pulsatile insulin release observed in vivo (28, 37). For this reason, an explanation for this type of oscillation in islets and β-cell lines is the main aim of this article.
It is critical to determine the degree of correlation between observed and simulated β-cell behavior. Several recent studies in the literature have shown additional experimental results consistent with the model that we present here, even though the authors' own interpretation is not identical to our model results. For example, KATP channel data (51) and data on glucose concentrations and oxygen consumption (2, 17) suggest that the ATP/ADP ratio may indeed oscillate on a time scale appropriate for the generation of slow oscillations. It was found that the ATP/ADP ratio drops rapidly when [Ca2+]i is raised and increases when [Ca2+]i falls (2, 19). These results are consistent with our calculations of [ATP]i behavior (Fig. 3). However, it should be pointed out that the [Ca2+]i-dependent changes in the ATP/ADP ratio are not driving oscillations in our model.
The assumption that [IP3]i can oscillate concomitantly with [Ca2+]i has received recent experimental support by the indirect demonstration of IP3 oscillations in ATP-stimulated Madin-Darby canine kidney epithelial cells (42). For murine β-cells in suspension, the intracellular IP3 level oscillated approximately in synchrony with changes in [Ca2+]i when driven by the sequential addition of glucose, an α2-adrenergic stimulus, and extracellular Ca2+ (5). It was also shown in RESULTS that our simulations are also consistent with data for [Na+]i oscillations.
All of these data are in reasonably good agreement with our model. However, it is important to emphasize that, although the modeled voltage and ion dependences are consistent with available data, the simultaneous measurements of the voltage dependence, relative conductance of different channels, and [Ca2+]i, [ATP]i, [Na+]i, [Ca2+]ER, and [IP3]i oscillations have not been directly demonstrated and need to be tested.
In addition, it should be pointed out that large oscillations in [ATP]i, [Ca2+]ER, [Na+]i, and [IP3]i should only be forthcoming in the slow [Ca2+]i oscillations, as was calculated above (Figs. 3, 5A, and 7). Fast [Ca2+]i oscillations can follow very small changes in other components (Figs. 4 and 5B), and it can be very difficult to measure these changes.
β-Cell Modeling: Identity of Slow Variables
Several types of models can explain electrical bursting and [Ca2+]i oscillations in β-cells. All known models have negative feedback from a slow variable driving the oscillations but differ in the identity of what the slow variable might be.
The first β-cell models were based on the hypothesis that the [Ca2+]i is the slow variable, influencing the membrane through a Ca2+-activated K+ current. However, the spectrophotometric measurements show that the time scale of the [Ca2+]i change is short relative to the oscillation period (see Refs. 26 and 78 for details). For this reason, we have not created a model in which [Ca2+]i is a possible driving slow variable.
The possibility that the ATP/ADP ratio is the slow variable driving electrical bursting by acting on the KATP channel was modeled by several investigators (46, 65, 82). However, the bursts and [Ca2+]i oscillations were principally obtained with a time constant appropriate for fast bursting in these models. Furthermore, it was shown in Refs. 74 and 75 that bursting and [Ca2+]i oscillations persist in the presence of a KATP channel blocker under conditions where IKATP is substantially decreased. This fact argues against a critical role for oscillations in KATP channels in these processes. We have also shown that, despite the fact that the change of mean ATP concentration can induce the oscillations or change their pattern, the ATP/ADP oscillation is not the driving force of [Ca2+]i oscillations in our model.
Several models have been proposed in which the dynamics of Ca2+ in the ER serves as the slow variable driving bursting (8, 13, 14, 26, 62). The ICRAN modifies the oscillation period and membrane potential in these models. We also followed these proposals for modeling fast oscillations. Our previous considerations (in RESULTS) showed that the dynamics of Ca2+ in the ER could be responsible for creating corresponding fast oscillations at a relatively high CaER2+ level, where there is a significant influence of [Ca2+]ER changes on ICRAN.
Although these models (8, 13, 14, 26, 62) have been applied mainly to an analysis of fast oscillations, this mechanism was also used for the generation of slow oscillations. Slow [Ca2+]i oscillations were obtained by Chay (14) at a decreased coefficient of Ca2+ER release from the ER. This retards the processes connected with CaER2+ dynamics and, in turn, transforms fast oscillations to slow. Increased [Ca2+]ER accompanies a decreased Ca2+ release from the ER. This means that the slow oscillations were generated at relatively high [Ca2+]ER in the model (14), contrary to our results. Emptying of the ER leads to termination of both fast and slow oscillations in the Chay model.
Arguments have been made against the CaER2+-dependent model for slow oscillations, since it seems to be at odds with data demonstrating that slow oscillations persist in the presence of Tg (6, 54, 64), an agent that blocks SERCA and empties the ER stores. These data were also analyzed above (in RESULTS), where we demonstrated that SERCA block can eliminate fast oscillations, while leaving intact the mechanism creating slow oscillations via the change in [Na+]i oscillations. For this reason, an expansion of the CaER2+-dependent model for fast oscillations to slow oscillations made by Chay (14) appears unjustified.
Recently, a Ca2+-sensitive K+ conductance was isolated that activates with a time constant of 2.3 s and inactivates with a time constant of 6.5 s in response to stimuli during voltage clamp (32). The authors concluded that that current, similar to a Ca2+-activated K+ current, played an important role in the generation of β-cell oscillatory activity. However, these time constants can only correspond to the fast oscillations.
Goforth at al. (31) have developed a new model to explain this slow K+ current, where [Ca2+]ER is also taken as the initial slow variable driving bursting. They introduced a new cytoplasmic Ca2+ pool localized between the ER and PM. Ca2+ released from the ER goes first to this intermediate Ca2+ pool. The changes of [Ca2+]i in this intermediate Ca2+ pool serve in this model to modulate the current through the Ca2+-activated K+ channels. We note also that the the value for the half-activation Ca2+ concentration assumed by Goforth at al. is much higher (0.7 μM) than that for the SK3 channels included here (0.1 μM).
However, Goforth at al. (31) chose this half-activation Ca2+ concentration arbitrarily to fit their model, specifically to address a potential high [Ca2+]i pool in the sub-PM region. Although such a Ca2+ pool has been proposed for excitatory cells, it has not been directly shown in β-cells, nor has the “real” value of the half-activation Ca2+ concentration for KCa channel in β-cells been determined. Kanno et al. (44) also suggest that this slow K+ current is determined not only by KCa channels. They found that activation and a reactivation of KATP channels can play a fundamental role in this slow K+ current in intact mouse pancreatic islets.
In addition, as was discussed above, every model of slow oscillations based on the proposal that [Ca2+]ER is the slow variable driving oscillations seems to be at odds with data demonstrating that slow oscillations persist in the presence of agents that block SERCA and empty the ER stores. This would restrict application of the Goforth et al. model to fast oscillations.
A very detailed model was developed by Magnus and Keizer (58, 59). They proposed (and simulated) that Ca2+ in mitochondria is the main regulator of the mitochondrial membrane potential, thus affecting ATP synthesis and corresponding [Ca2+]i oscillations. However, the view that mitochondrial change of calcium regulates ATP production under physiological conditions has been questioned (2, 27), and experimental support is still incomplete for the idea that changes in the ATP/ADP ratio are the main driving mechanism for all Ca2+ oscillations (see above).
Despite the clear evidence that IP3 is one of the critical cellular messengers inducing Ca2+ mobilization from intracellular stores (38), IP3 kinetics have not previously been included in the models of the pancreatic β-cell. Our model produced a simulation with a good match to data reported from several laboratories. The simulations also underlined how the IP3 concentration was critical to understanding the pattern of slow [Ca2+]i oscillations even though IP3 does not account for the mechanism of oscillations.
Our simulations suggest that the change in [Na+]i is one of the most important slow variables that govern slow [Ca2+]i oscillations. Kitasato et al. (47) proposed such a role for Na+. According to these authors, an increase in [Ca2+]i results in an increase in [Na+]i because of Na+/Ca2+ exchanger activity. The increase in cytosolic Na+ results in enhanced ATP consumption at Na+/K+ pump sites and in a decreased ATP/ADP ratio, which in turn could lead to activation of KATP channels followed by membrane repolarization. A decrease in [Ca2+]i results in the opposite order of events during a silent phase. This scenario was modeled by Miwa and Imai (65) as the underlying mechanism of [Ca2+]i oscillations. However, this mechanism does not explain the experimental data suggesting that ATP/ADP ratio changes are not the main driving force for [Ca2+]i oscillations (see above).
Our model simulates an increase in [Na+]i and a corresponding decrease in [ATP]i during the depolarized phase by use of a similar mechanism. However, our model shows that an increase of [Na+]i during the depolarized phase and a corresponding increase of outward current through Na+/K+ pumps could be sufficient for PM repolarization that leads to transition to silent phase, even though a change of ATP concentration is not directly driving oscillations.
Slow oscillations of PM potential and [Ca2+]i were obtained in previous models (e.g., Refs. 7 and 78) by introducing a special slow variable that does not have a clear metabolic basis. Our work is an attempt to find some real metabolic basis for this suggestion. According to our analysis, the Na+ concentration in the cytoplasm may be the desired slow variable.
We also found that some coefficients used in previous models do not always correspond to physiological mechanisms. It has been suggested in previous models that Ca2+ is removed from the cytoplasm only by sequestration mechanisms that do not create any PM currents. For example, in one model (13), the mitochondrion sequesters intracellular Ca2+ with a rate constant of 3,000 s-1, and this calcium is not recycled to the cytoplasm. In other examples, the Ca2+ removal rate through the PM (e.g., by sequestration in secretory granules) was assigned the following coefficients: 120 s-1 (82), 190 s-1 (58), and 640 s-1 by Gall and Susa (26) in their models II and III. This means that, in all of these models, total cytoplasmic Ca2+ disappears from the β-cells in only seconds as the result of mechanisms that do not create any current through the PM. We could not find any data to support such a significant flux created by sequestration mechanisms. Ca2+ is removed mainly by well-known PM Ca2+ pumps in our model that create a corresponding PM current (Eq. 8). The coefficient for the sequestration rate of Ca2+ through the PM (ksg in Eq. 24) was estimated in our models as 0.1 s-1, and we consider this to be a more reasonable value.
In the model by Miwa and Imai (65), only Na+/K+ active transporters and the Ca2+ pumps on the PM were employed to simulate ATP consumption in β-cells. Furthermore, the rate constant of ATP production from ADP was ∼0.03 s-1 at high glucose levels. Other ATP consumption mechanisms exist in cells, which we simulate in our model (Eq. 27). The rate constant of ATP production from ADP was also considerably increased in our model (up to 10 times in Fig. 3) to be in agreement with experimental observations of ATP production rate (see MATERIALS AND METHODS). Furthermore, no flux analysis or analysis of oscillations at fixed values of ATP concentrations was made in (65). On the basis of simulations in our model, it seems possible that, in the model of Miwa and Imai, a change of Na+ concentration also can directly effect PM potential by a change of INa/K.
Several additional slow variables can also take part in the regulation of oscillations in β-cells. For example, the concentrations of such ions as K+ and Cl- were used as independent variables in the model (65). Unquestionably, K+ concentration in the cytoplasm is very important. However, all cells, including β-cells, have 100–130 mM K+ concentration in the cytoplasm, and its relative change must be small, particularly during fast oscillations (65). For this reason, we did not include K+ cytoplasmic concentration as an independent variable in this model. It has also been found that glucose stimulation of the β-cells is coupled to an increase in their Cl- permeability and that the oscillatory Ca2+ signaling may be critically dependent on transmembrane Cl- fluxes (22). It seems plausible that additional kinetic equations for K+,Cl-, and other ions and their channels may become necessary for this model in the future to include emerging aspects of β-cell physiology.
In conclusion, our results support the hypothesis that it is necessary to include all essential ionic channels and pumps that have been identified in β-cells to employ mathematical models to correctly simulate the set of experimental results. Results from the simulations support the hypothesis that [Ca2+]ER and [Na+]i can be the main variables driving both fast and slow [Ca2+]i oscillations and that the regulation of Ca2+ leak from the ER by IP3 is also important to explain the pattern of slow calcium oscillations. The results of the modeling favor a scenario where two or more slow variables interact to produce oscillations, in agreement with other recent models (6, 7). On the basis of our analysis, it seems likely that differences in relative activity of channels and pumps could explain the variations in the electrical bursts and [Ca2+]i oscillations between the islets, cell lines, and single β-cells.
APPENDIX A: NA+/K+ ACTIVE TRANSPORT
The general model (12) was employed in a modified form according to Miwa and Imai (65). A1 where where PNaK is a coefficient for resulting current in the presence of saturating level of ATP; [P] is an inorganic phosphate concentration; b1, b2, b3, b4, b*5, b6, f1, f2, f3, f*5, and f6 are the rate coefficients: [P] = 4,950 μM, f1 = 2.5 10-10 μM-3 ms-1, f2 = 10 ms-1, f3 = 0.172 ms-1, f4 = 1.5 10-8 μM-2 ms-1, f*5 = 0.002 μM-1 ms-1, f6 = 11.5 ms-1, b1 = 100 ms-1, b2 = 0.0001 μM-1 ms-1, b3 = 1.72 10-17 μM-3 ms-1, b4 = 0.0002 μM-1 ms-1, b*5 = 0.03 ms-1, b6 = 6 10-7 μM-1 ms-1).
APPENDIX B: GENERAL GATING VARIABLE FOR IKDR
According to Gall and Susa (26) and Sherman (78), the general gating variable for IKDr can be determined as B1 where B2 B3 where n∞ is the steady-state value of the general gating variable, Vn is the half-activation potential, Sn is the slopes at half-maximal potential, τn(V) is the relaxation time constant, a, b, c, and Vτ are the adjustable coefficients. The coefficients were fitted on the basis of previous models (26, 82): Vn = -14 mV, Sn = 7 mV, a = 65 mV, b = 20 mV, c = 20 ms, and Vτ = -75 mV.
APPENDIX C: FRACTION OF KATP CHANNELS OPEN
For the fraction of KATP channels open, we adopt a detailed kinetic model that was created by Magnus and Keizer (58) and that was based on the results of Ref. 43. The equation for the fraction of open channels can be included on the basis of that model if we take into account the parameters determining adenine nucleotide concentrations (see Table 2 from Ref. 58). Then the resulting fraction of KATP channels open has the form C1 OKATP is the fraction of channels open. Kdd, Ktd, and Ktt are the dissociation constants (Kdd = 17 μM, Ktd = 26 μM, and Ktt = 1 μM). According to Hopkins et al. (43), Eqs. 19 and C1 reproduce well the experimental data.
We thank A. Kuznetsov and M. Roe for helpful discussions. This work was partially supported by National Institute of Diabetes and Digestive and Kidney Diseases Grants DK-44840, DK-20595, and DK-48494, and the Blum-Kovler Foundation.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2003 by American Physiological Society