|Home | About | Journals | Submit | Contact Us | Français|
Pancreatic beta-cells respond to rising blood glucose by increasing oxidative metabolism, leading to an increased ATP/ADP ratio in the cytoplasm. This leads to a closure of KATP channels, depolarization of the plasma membrane, influx of calcium and the eventual secretion of insulin. Such mechanism suggests that beta-cell metabolism should have a functional regulation specific to secretion, as opposed to coupling to contraction. The goal of this work is to uncover contributions of the cytoplasmic and mitochondrial processes in this secretory coupling mechanism using mathematical modeling in a systems biology approach.
We describe a mathematical model of beta-cell sensitivity to glucose. The cytoplasmic part of the model includes equations describing glucokinase, glycolysis, pyruvate reduction, NADH and ATP production and consumption. The mitochondrial part begins with production of NADH, which is regulated by pyruvate dehydrogenase. NADH is used in the electron transport chain to establish a proton motive force, driving the F1F0 ATPase. Redox shuttles and mitochondrial Ca2+ handling were also modeled.
The model correctly predicts changes in the ATP/ADP ratio, Ca2+ and other metabolic parameters in response to changes in substrate delivery at steady-state and during cytoplasmic Ca2+ oscillations. Our analysis of the model simulations suggests that the mitochondrial membrane potential should be relatively lower in beta cells compared with other cell types to permit precise mitochondrial regulation of the cytoplasmic ATP/ADP ratio. This key difference may follow from a relative reduction in respiratory activity. The model demonstrates how activity of lactate dehydrogenase, uncoupling proteins and the redox shuttles can regulate beta-cell function in concert; that independent oscillations of cytoplasmic Ca2+ can lead to slow coupled metabolic oscillations; and that the relatively low production rate of reactive oxygen species in beta-cells under physiological conditions is a consequence of the relatively decreased mitochondrial membrane potential.
This comprehensive model predicts a special role for mitochondrial control mechanisms in insulin secretion and ROS generation in the beta cell. The model can be used for testing and generating control hypotheses and will help to provide a more complete understanding of beta-cell glucose-sensing central to the physiology and pathology of pancreatic β-cells.
The appropriate secretion of insulin from pancreatic β-cells is critically important for energy homeostasis. Pancreatic β-cells are adapted to sense blood glucose and other secretagogues to adjust insulin secretion according to the needs of the organism. Rather than activating specific receptor molecules, glucose is metabolized to generate downstream signals that stimulate insulin secretion. Pancreatic β-cells respond to rising blood glucose by increasing oxidative metabolism, leading to increased ATP production in mitochondria and in an enhanced ratio of ATP to ADP (ATP/ADP) in the cytoplasm [1-3]. The increase in intracellular ATP/ADP closes the ATP-sensitive K+ channels (KATP), decreasing the hyperpolarizing outward K+ flux. This results in depolarization of the plasma membrane, influx of extracellular Ca2+ through the voltage-gated Ca2+ channels, a sharp increase in intracellular Ca2+ and activation of protein motors and kinases, which then mediate exocytosis of insulin-containing vesicles [2-5]. The currently accepted processes of glucose metabolism and Ca2+ handling in the cytoplasm and mitochondria of β-cells considered in this analysis are summarized in Figure Figure11[1-4].
A brief summary of these processes includes the following steps. Glucose enters β-cells by facilitated diffusion through glucose transporters (GLUT1 and 2). While this process is not limiting in β-cells , the next irreversible step, glucose phosphorylation, is catalyzed by a single enzyme, glucokinase (GK). This enzyme is specific for metabolic control in the β-cell and hepatocyte, because the Km of GK for glucose is ~8 mM, a value that is almost two orders of magnitude higher than that of any other hexokinase. This step appears to be rate limiting for β-cell glycolytic flux under normal physiological conditions, so that GK is regarded as the β-cell 'glucose sensor' [1,3], underlying the dependence of the β-cell insulin secretory response to glucose in the physiological range.
Pyruvate is the main end product of glycolysis in β-cells and essential for mitochondrial ATP synthesis. In the mitochondrial matrix, pyruvate is oxidized by pyruvate dehydrogenase to form acetyl-coenzyme A (acetyl-CoA). Acetyl-CoA enters the tricarboxylic acid (TCA) cycle to undergo additional oxidation steps generating CO2 and the reducing equivalents, flavin adenine dinucleotide (FADH2) and NADH. Oxidation of reducing equivalents by the respiratory chain is coupled to the extrusion of protons from the matrix to the outside of the mitochondria, thereby establishing the electrochemical gradient across the inner mitochondrial membrane (Figure (Figure1).1). The final electron acceptor of these reactions is molecular oxygen, as in other eukaryotic cells. The electrochemical gradient then drives ATP synthesis at the F1F0-ATPase complex to phosphorylate mitochondrial ADP, thereby linking respiration to the synthesis of ATP from ADP and inorganic phosphate (Figure (Figure1).1). Adenine nucleotide translocase (ANT) exchanges matrix ATP for ADP to provide ATP for energy consuming processes in the cytosol. Some cytosolic ATP is also produced in the latter part of glycolysis. However, this appears to be of minor consequence relative to that subsequently generated in the mitochondria, which represents an estimated 90% of the total β-cell ATP production [7,8].
The cytoplasmic Ca2+ signal is coupled to mitochondrial Ca2+ handling (Figure (Figure1).1). The balance of Ca2+ influx and efflux determines the matrix Ca2+ level involving the Ca2+ uniporter and the mitochondrial Na+/Ca2+ exchanger, respectively. Ca2+ influx into mitochondria is amplified by hyperpolarization of the inner mitochondrial membrane [9,10]. Inside the organelle, Ca2+ activates several matrix dehydrogenases (for example, pyruvate dehydrogenase). Mitochondrial Ca2+ may also directly stimulate ATP synthase . The nutrient-dependent Ca2+ rise in the cytosol further activates ATP hydrolysis [7,10,12,13].
An important β-cell specialization is the very low expression of lactate dehydrogenase (LDH), the enzyme catalyzing the conversion of pyruvate to lactate [1,14,15]. A low level of LDH expression in insulin-secreting cells is important to preferentially channel pyruvate towards mitochondrial metabolism (see [1,10,16]). However, the low LDH levels likely leads to activation of compensatory mechanisms because NAD+-dependent glycolytic enzymes (e.g., glyceraldehyde 3-phosphate dehydrogenase) require that cytoplasmic NADH must be re-oxidized to NAD+. This reaction is usually catalyzed by LDH, but because β-cells cannot use this pathway effectively, these cells must re-oxidize cytoplasmic NADH by activation of two mitochondrial hydrogen shuttles (Figure (Figure1),1), the malate-aspartate shuttle and the glycerol phosphate shuttle [15,17-19].
Glucose signaling in β-cells has several other peculiarities, including generation of multiple oscillations in metabolism, mitochondrial membrane potential Ψm) and NADH, mitochondrial and cytoplasmic Ca2+ and, ultimately, the oscillations of insulin secretion [5,20-23]. The coupling of these various oscillators is not clearly understood. In addition, the respiratory rate is lower and relative leak activity is higher in isolated β-cell mitochondria (as found in a cultured β-cell line) compared with isolated mitochondria from skeletal muscle [24,25]. These observations need clarification to better understand how mitochondrial processes are linked with insulin secretion.
The unique character of the β-cell response to glucose is usually attributed solely to glucokinase. Because of its near-dominant control of glycolytic flux, this enzyme is thought to govern the ATP/ADP ratio and insulin secretion almost exclusively [1,3]. While glucokinase certainly exerts a critical level of control on downstream events, other cytoplasmic and mitochondrial processes also play an essential role in glucose-stimulated insulin secretion (GSIS) [1,2,10]. In particular the relatively high flexibility of the ATP/ADP ratio in β-cells may be accounted for, at least partly, by mitochondrial peculiarities as well as by properties of glucokinase [24,26,27]. For these reasons it is critical to develop a comprehensive understanding as to how cytoplasmic and intramitochondrial fuel metabolism is coupled to fuel availability and thereby "sensed."
The goal of this work is to determine the contribution of the cytoplasmic and mitochondrial processes regulating GSIS using a mathematical modeling approach. Mathematical modeling can be a powerful systems biology tool allowing quantitative descriptions of the control individual components exert over the whole biological system. Several mathematical approaches in the literature have provided quantitative estimates of energetic and mitochondrial processes in pancreatic β-cells. However, these models are limited in the pathways that are considered, so that a more comprehensive approach is now necessary.
The first detailed β-cell model was developed by Magnus and Keizer [28-30]. However, several mechanisms used for simulations in this model have recently been reevaluated. For example, steady-state electron transport and the F1F0 ATPase proton pump were modeled according to the "six states proton pump mechanism" . This mechanism does not correspond to the present understanding of the function of the electron transport chain (ETC) and the mitochondrial F1F0 adenosine trisphosphatase (see for example ). Models of the LDH and NADH shuttles were not included, and mitochondrial fluxes may also have been overestimated in this model (see below). The main goal of these models was to examine the possible mechanisms underlying oscillations in pancreatic β-cells, not biochemical regulation of β-cell glucose sensitivity that we are focused on here.
A complex kinetic model of the metabolic processes in pancreatic β-cells based on in vitro enzyme kinetics was recently developed . However, while heroically complicated models with numerous parameters and enzyme activities are interesting, they require data on in vivo enzyme activities and coefficients that are not readily available.
Enzyme activity measurements in vitro, often used in models, may not reflect enzyme activity in vivo . For example, experimental kinetic data for isolated mitochondria and the parameters evaluated for mitochondrial processes from experiments with intact cells may differ significantly . For these and other reasons previous models of pancreatic β-cell energetics and mitochondrial calcium regulation fall short of a comprehensive explanation of the mechanisms of β-cell sensitivity.
To address this we have developed a specific quantitative, kinetic model (see Appendix) of the core processes of β-cell cytoplasmic and mitochondrial energetic based on a simplified map of the biochemical pathways schematized in Figure Figure1.1. We included the most recent experimental characterizations of the majority of processes in the model to insure accuracy. However, for simplification, we modeled only those regulatory couplings that we have deemed most crucial for the β-cell metabolic regulation based on experimental evidence. The model includes the dynamic equations for cytoplasmic ADP, NADH and glyceraldehyde 3-phosphate, mitochondrial Ψm and NADH, mitochondrial and cytoplasmic Ca2+ and pyruvate. When available we used the values of the coefficients determined for living cells rather than for isolated enzymes and cell-free mitochondria (Appendix).
We show that this model has qualitative properties consistent with expectations for the pancreatic β-cell including showing appropriate oscillations in mitochondrial metabolism and Ca2+ concentration. The model also reproduces simultaneous measurements of the behavior of multiple constituents within the cytoplasm and mitochondria such as NADH, Ca2+ and Ψm at high temporal resolution. We also discuss specific differences in muscle and β-cell mitochondrial function, providing insight into essential control properties of the β-cell. Furthermore, predictions on the dynamics of as yet unmeasured molecules could be made, and the model further tested by verifying these predictions.
Nutrient-stimulated insulin secretion in β-cells is impaired in the diabetic state. This may result from impaired glucose-induced ATP/ADP ratio elevation in β-cells [26,35]. Furthermore, it is becoming increasingly clear that the development of type 2 diabetes is associated with mitochondrial dysfunction [27,35-37]. Insulin signaling also effects mitochondrial function in β-cells . Thus, knowledge of the mechanisms of regulation of ATP production and consumption are central to understand β-cell glucose-sensing and mechanisms of dysfunction in type 2 diabetes.
The model was used to simulate data obtained using several experimental protocols. Under resting low glucose concentration the simulated values of the cytoplasmic and mitochondrial variables are consistent with experimentally reported data as indicated in Table Table1.1. Then the model was used to examine the steady-state changes of the state variables and fluxes. Figure Figure22 shows the responses of model simulation to steps in the glucose concentration observed in successive steady-states. Glucokinase catalyzed the rate-limiting step of glycolysis with a steep dependence on glucose concentration in the range 4-25 mM. Enhancement of glucose concentration led to an increase in glycolytic flux, glyceraldehyde 3-phosphate (G3P) and pyruvate concentrations (Figure 2A,B). This process accelerated pyruvate reduction and decarboxylation leading to increased [NADH]m (Figure (Figure2B).2B). [NADH]m was oxidized by the ETC, raising the rate of mitochondrial O2 consumption. Oxidation of mitochondrial NADH by the respiratory chain increased the membrane potential directly via protons pumped out of the matrix. Ψm was dissipated by proton-leak reactions and the activity of the phosphorylation apparatus, which included the phosphate carrier and the ATP synthase (Figure (Figure1).1). The net result of these processes was establishment of an elevated Ψm. The hyperpolarization of the inner mitochondrial membrane resulted in increased ATP production by F1F0 ATPase, decreased [ADP]c and a corresponding increased ATP/ADP ratio (Figure (Figure2D).2D). The phosphorylation rate (Jph) reached saturation at high glucose concentration (Figure (Figure2C)2C) as a consequence of decreased [ADP]c and saturated Ψm(see Equation 7 and Figure Figure1212 in Appendix). Simultaneously, [Ca2+]c increased with increased ATP/ADP ratio according to the empirical Equation 23 (Appendix). We also simulated the steady-state response of free mitochondrial matrix Ca2+ to changes in cytoplasmic Ca2+ concentration, Ψm and finally glucose (Figure (Figure2E2E).
As expected, our simulations were consistent with experimental data. Glucose utilization increased lactate synthesis, O2 consumption and CO2 production [1,7,39,40]. Both cellular [G3P] and [PYR] increased after simulating increased extracellular glucose (Figure (Figure2B).2B). This result is consistent with the finding of increased glycolytic intermediates and pyruvate after glucose challenge in the INS1 β-cell line  as well as the increase in Ψm with increased glucose in mouse islets [20-22,42-44].
An accurate measurement of lactate output in β-cells from isolated islets is difficult to obtain because LDH expression in non-β-cells is considerably higher than in β-cells, and high rates of lactate output may also originate from cells in the centers of isolated islets that are prone to oxygen depletion and necrosis [39,45]. However, the oxidative production of CO2 from [3,4-14C]glucose represented close to 100% of the total glucose utilization in purified rat β-cells  indicating that lactate output should not exceed several percent. Very low lactate output was also found in β-cell lines . Our simulated small lactate output in Figure Figure2A2A is consistent with these experimental data.
The results of the simulation (Table (Table11 and Figure Figure2B)2B) were also consistent with the range of measured [NADH]m reported previously. For example, the concentration of free NADH in mitochondria of intact pancreatic islets at resting glucose levels (4-5 mM) is about 60 μM and the maximum mitochondrial glucose-induced increase in free NAD(P)H reached 75 μM . The simulated increased [Ca2+]c versus glucose concentration (Figure (Figure2E)2E) was also in agreement with previous reports (see for example [42,48-50]).
Several studies have confirmed an increase in the ATP/ADP ratio in response to high glucose (see for e.g. [3,7,12,26,51]. A simultaneous rise in ATP/ADP and NADH/(NAD+ + NADH) ratio was found in rat islets , and NAD+/NADH was increased in rat β-cells and in the MIN6 β-cell line in response to high glucose . The rise in ATP/ADP ratio as well as in relative NAD(P)H, Ψm, [Ca2+]m and oxygen consumption were also observed with glucose stimulation in control INS-1 cells [54,55]. Our simulations are generally consistent with these data.
The model suggests a possible reconciliation of several apparent contradictions between live cell experimental data and regulation of mitochondrial energetics obtained in experiments with isolated mitochondria.
1. A basic principle of mitochondrial energetics is given by the inverse relationship between the respiratory flux and Ψm, i.e. the higher Ψm, the lower the respiration rate [24,56]. We simulated this relationship using Equation 5C (Appendix). However, the electron transport rate (Jhres) and O2 consumption increased simultaneously with Ψm in our simulation of the β-cell (see Figure Figure2C)2C) as well as in vivo (see above). The model offers the following explanation of this contradiction. In our model (as in living cells) the electron transport rate (Equation 5) depends on at least two factors: one is a decrease in the electron transport rate with an increase in Ψm (Equation 5C) but another factor is an increase of this rate with increased substrate concentration (NADH) (Equation 5A). Increasing the electron transport rate simultaneously with Ψm means that the enhancement of Jhres as a result of the increased [NAPH]m was greater than its decrease with the rise of Ψm following the step increase in glucose. Substrate concentrations are usually maintained at constant or saturated levels in experiments with isolated mitochondria, where one can only see inhibition of the electron transport rate with increased Ψm.
2. The respiratory control hypothesis for ATP production in intracellular mitochondria was based on experiments with isolated mitochondria which found that ADP availability to the ATP-synthase is the limiting factor for mitochondrial ATP production , that is, the rate of ATP synthesis should decrease with decreased [ADP]c. This mechanism corresponds to Equation 7A in our model. Experimentally this hypotheses has been tested in permeabilized clonal β-cells, where ATP/ADP ratios can be externally fixed showing that a decrease in [ADP] led to decreased O2 consumption . However, an increased ATP/ADP ratio (usually due to decreased [ADP]c) coincidentally with increased respiration rate and oxidative phosphorylation has been firmly established for pancreatic β-cells as a signal for GSIS in response to increased glucose [1,3,7,12,26,51,55]. Similar results were obtained in our simulation of β-cell shown in Figure Figure2D.2D. At first glance these data seem inconsistent with the expected inhibition of respiration with decreased ADP concentration [3,55].
Our analysis resolves this apparent contradiction. In our model the ATP synthesis rate is dependent on at least two factors: one is a decreased ATP synthesis rate with decreased [ADP]c (Equation 7A) but another factor is an increased ATP synthesis rate with increased Ψm(Equation 7B). Our simulation shows that enhancement of ATP production with increasing Ψm was greater than its decrease as a result of decreasing [ADP]c following a step increase in glucose. As a result, ATP synthesis and respiration rate increase despite decreased [ADP]c and the ATP/ADP ratio increased with a step glucose increase (Figure (Figure2D).2D). These simulations imply that glucose challenge can lead to simultaneous increases in Ψm, the ATP/ADP ratio and in the rates of mitochondrial ATP synthesis and respiration.
The concentrations of free ATP and ADP in the cytoplasm were used in our model since only free molecules can take part in reactions. However, the free ATP concentration is close to its total concentrations, whereas the fraction of bound ADP may be substantial [58,59]. On the other hand, most estimated ATP/ADP ratios are based on measurements of total nucleotide content [7,12,51]. For this reason, the measured ATP/ADP ratio of total ATP and ADP nucleotide content is likely to be substantially smaller than ratio of concentration of the free components, simply because the measured total ADP content includes bound ADP. Therefore, it is not surprising that the simulated ATP/ADP ratio change in Figure Figure2D2D using free nucleotide concentrations is greater than that in published experimental data (see for example [7,12,51].
According to our simulation only a small increase in the ATP concentration occurred following glucose challenge (not shown). A decrease in the free [ADP]c is the main factor leading to an increase in the ATP/ADP ratio following increased glucose in Figure Figure2D.2D. This simulation is in agreement with experimental data and can be a consequence of the initial high ATP/ADP ratio even with a low glucose level in our model (see Table Table1).1). For this reason, the ATP concentration cannot be increased significantly if the total adenine nucleotide concentration is kept constant, whereas the relative [ADP] may undergo a pronounced decrease (see our previous publication  for a detailed consideration of this question).
β-cell regulatory mechanisms endow this cell type with unique metabolic properties to control insulin secretion in comparison with metabolism in other cell types. For example, liver cells maintain a stable ATP/ADP equilibrium while respiring at widely varying rates . Cardiac myocytes can increase, by three- to sixfold, the rate of cardiac power generation, myocardial oxygen consumption, and ATP turnover in the transition from rest to intense exercise . Nevertheless, at high work states the myocardial ATP and ADP concentrations are maintained at a relatively constant level despite the increased turnover rates [34,62].
Specific β-cell respiratory mechanisms can be illustrated by comparing isolated mitochondria from skeletal muscle and cultured β-cells. The rate of respiration was higher (>5.5 fold) and the relative leak rate was significantly lower at any Ψm value in isolated mitochondria from skeletal muscle than in those from cultured β-cells [24,25]. We examined how these differences effect mitochondrial function by simulating the conditions of work in muscle mitochondria (Figure (Figure3).3). Mitochondrial NADH and cytoplasmic ADP concentration are maintained at a relatively high and constant level in muscle cells [34,62]. To simulate this, the concentration of [NADH]m was set as a constant reflecting this concentration for high glucose level in a β-cell (25 mM). [ADP]c was also set to an elevated constant level (700 mM), that was 5-fold higher than the calculated [ADP]c level (at 9 mM glucose) in β-cells.
Figure Figure33 shows the results of simulations in which the maximal rate of ETC (Vme) was increased in steps. Mitochondrial F1F0 ATPase activity (Vmph) was unchanged. Simulated Ψm and the rate of ATP production (Jph) were significantly increased with an increased Vme, such that F1F0 ATPases work with maximal activity under these conditions (compare Jph in Figure Figure2C2C and Figure Figure3).3). This can be explained by the high Ψm (more electronegative) as well as by the increased [ADP]c in simulated muscle cells in comparison with β-cells. Note that the rate of ATP production (Jph) depended only slightly on Ψm change when Ψm was increased above 160 mV, since these levels of Ψm were saturating for F1F0 ATPase activity (Appendix, Figure Figure1212).
This indicates that the F1F0 ATPase can work in muscle cells with maximal productivity during increased respiration activity because ADP concentration and Ψm are supported at relatively higher levels. It thus appears that a decrease in the efficiency of mitochondrial energy production with decreased Ψm can lead to a relatively high degree of control on the phosphorylation potential in β-cells, i.e. a change in Ψm leads to a large change in Jph. Interestingly, the simulated relative leak (Jh1) magnitude was significantly lower in the muscle cell simulation in comparison with respiration rate (evaluated as Jhres) at increased Vme even with invariant coefficients for the proton leak, since the rates of respiration and ATP production were highly increased but a coefficient of leak (Jh1) would remain as constant (Figure (Figure33).
Our simulations help explain the data of Affourit and Brand [24,25] showing decreased respiratory and increased relative leak activity in isolated β-cell mitochondria. This suggests that mitochondrial glucose sensitivity in β-cells results from decreased respiratory activity compared with F1F0 ATPase activity. This leads to mitochondrial work at decreased Ψm that is in the region where variations in Ψm should result in an increased sensitivity to glucose. Decreased respiratory activity in β-cells leads to a decreased ATP production rate by the F1F0 ATPase. However, this gives β-cells the ability to adaptively change the ATP/ADP ratio in response to changes in glucose concentration.
Interestingly, the oxidative phosphorylation rate (per g dry weight) was significantly lower in pancreatic islets even in high glucose, compared with brain or heart , supporting our suggestion regarding decreased respiratory activity in β-cells.
Our model features increasing proton-leak with increased Ψm (Equation 8, Appendix). We simulated how changes in leak activity affect the response of the variables. To do this we increased the regulated coefficient of proton leak (P1r in Equation 8, Appendix) threefold leading to an increase in the total proton leak rate by twofold (Figure (Figure4).4). As expected, one effect was to reduce the inner membrane potential that causes a corresponding right-shift in the ATP/ADP ratio and [Ca2+]c response to glucose. To simulate decreased leak activity (for example following decreased uncoupling protein expression) we set the regulated proton leak coefficient equal to zero in Equation 8 (P1r = 0) (Figure (Figure4).4). In this way the general leak activity was decreased by 50%. Decreased leak activity increased Ψm and reduced the sensitivity range of the inner membrane potential to glucose leading to a left-shift in the ATP/ADP ratio and [Ca2+]c response.
These simulations show that proton leak can modulate GSIS by shifting the dependence on glucose of the ATP/ADP ratio and [Ca2+]c, altering cellular sensitivity to glucose challenge. This effect of proton leak is only possible when the ATP/ADP ratio can be regulated by changes in Ψm, i.e. when Ψm lies below the β-cell maximal level for ATP production. On the other hand, in muscle cells Ψm can be maintained at a high level (see above) and changes in Ψm exert an insignificant effect on the ATP production rate.
This role of uncoupling agents can be illustrated by considering the experimental data for the principal β-cell uncoupling protein 2 (UCP2) (see for review [63,64]. For example, overexpression of UCP2 in normal rat islets diminished the change in mitochondrial membrane potential in response to glucose, reduced cytoplasmic ATP content, GSIS and mitochondrial ROS production [65,66]. β-cells exposed to free fatty acids displayed a lower mitochondrial membrane potential (less electronegative) and a decreased glucose-induced hyperpolarization. These effects were due to increased activity of UCP2 . Conversely, UCP2-deficient mice demonstrated increased ATP production and improved GSIS . In pancreatic islets from wild type but not Ucp2-knockout mice, genipin, a cell-permeant compound that was reported to inhibit UCP2-mediated proton leak, increased the mitochondrial membrane potential and cytosolic ATP, closed KATP channels, and stimulated insulin secretion .
Our simulation gave similar results (Figure (Figure4),4), showing that at a constant glucose level, increased uncoupling protein activity leads to decreases in Ψm, the ATP/ADP ratio and [Ca2+]c. However, our simulation also shows a novel aspect of this problem: an ability of uncoupling proteins (and other uncoupling agents) to shift the glucose dependence of the ATP/ADP ratio and [Ca2+]c. This shifting mechanism is not simply varying the rates of a few processes. In our model, the shifting mechanism is generally characterized by small changes in Ψm, ATP/ADP ratio or [Ca2+]c at low and high glucose levels. However, larger changes of [Ca2+]c can be expected in the region of physiological glucose concentrations (Figure (Figure4),4), a hypothesis that needs to be further tested.
The shifting set-point mechanism would regulate insulin secretion particularly during a fluctuating nutrient supply. For example, UCP2 has been shown to be both induced and activated by exposure of rodent islets to high free fatty acids that cause mild uncoupling [64,67,69,70]. According to our analysis uncoupling can lead to a right shift in the glucose dependence of the ATP/ADP ratio and [Ca2+]c. This shift can result in decreased insulin secretion at moderate levels of glucose and at high levels of free fatty acids in blood. It may be a mechanism for restricting glucose consumption under conditions of increased serum fatty acid concentrations, when muscle cells can use free fatty acids as fuel rather than glucose.
Interestingly, this may be a physiologically important mechanism to blunt insulin production in starved animals even without increasing free fatty acid concentration. During starvation homeostatic mechanisms attempt to maintain a minimal acceptable blood glucose concentration, in part to maintain neurons that cannot use free fatty acids. Starvation induces UCP2 expression and reduces cellular NADH generation in response to glucose in mouse β-cells that would limit insulin secretion to reduce glucose uptake in muscle and adipose cells .
Our analysis supports an important role for uncoupling agents in β-cells that can coordinate the appropriate response of β-cells to fluctuating nutrient supply (see for example [24,63]. However, a pathological side effect might also accompany shifting ATP/ADP ratio and [Ca2+]c sensitivity to glucose, because decreased sensitivity of insulin secretion in response to glucose is a characteristic property of Type 2 diabetes [27,37]. Simulation with this model supports the idea that either an under or over-expression of UCP2 may lead to a failure of β-cells to properly respond to glucose.
Since β-cells have low levels of LDH (see Introduction) [1,14,15], we included a low level of LDH activity under basal conditions (Appendix). This led to a small lactate flux within the cytoplasm at low and medium glucose levels in our simulations (Figure (Figure2A).2A). However, the simulated lactate production increased significantly in response to high glucose. This can be accounted for by increased [NADH]c/[NAD+]c (see Equation 3, Appendix).
We also simulated a rise in LDH activity (Figure (Figure5).5). This led to accelerated conversion of pyruvate to lactate, that decreased NADH production in mitochondria from pyruvate and the corresponding [NADH]m, Ψm, ATP/ADP ratio as well as [Ca2+]c in response to increased glucose. This was all a consequence of the increased fraction of the glycolytic flux that was directed to lactate production.
Our simulations confirm that low levels of LDH expression in insulin-secreting cells are important for the correct channeling of pyruvate towards mitochondrial metabolism (see [1,16]). However, we also found that net lactate production increases significantly when extracellular glucose is increased (see Figures Figures22 and and5).5). For this reason, even low LDH activity can be an effective safeguard to prevent mitochondrial overexcitation at high glucose levels, where [Ca2+]c concentration is already saturated and increased Ψm can lead to increased ROS production (see below).
Interestingly, overexpression of LDH in single MIN6 β-cells diminished their response to glucose as measured by mitochondrial NAD(P)H, Ψm, cytosolic free ATP and Ca2+ and led to a right shift in the glucose response of insulin secretion , all of which are simulated by our model. Since islet levels of LDH were increased in a rat pancreatectomy model of type 2 diabetes  and β-cell lines such as INS-1 have increased LDH activity that can partially explain their decreased sensitivity to glucose , overexpression of LDH may also be a possible mechanism of β-cell failure in specific cases.
The free cytoplasmic [NAD+]c/[NADH]c ratio can vary from about 700:1 to 200:1 in rat liver . Zhang et al.  measured free [NAD+]c/[NADH]c in Cos7 cells and found a ratio of 644. Patterson et al.  found that the NAD(P)H concentration (NADH plus NADPH) was much lower in the cytoplasm than in mitochondria in pancreatic islets. Our calculations for the basal [NAD+]c/[NADH]c ratio level (Figure (Figure2B)2B) yield a magnitude of about 750:1 (at a glucose of 9 mM). Generally, our simulated [NAD+]c/[NADH]c ratios were in the range of experimentally observed values. As free NAD+ levels greatly exceed those of NADH, large changes in the [NAD+]c/[NADH]c ratio do not require correspondingly large changes in free [NAD+]c. Thus, changes in cytoplasmic redox level could be manifested primarily through variation in [NADH]c.
Pyridine nucleotides are not only key molecules for metabolic conversions, but also can serve as critical signaling molecules, such as the cytoplasmic NAD(P)H/NAD(P)+ ratio [2,42,53,75]. Here we have focused only the energetic aspects of pyridine nucleotides and modeling. Generally, redox shuttles are responsible for maintaining and restoring cytoplasmic NAD(P)H/NAD(P) + ratios and the cytoplasmic/mitochondrial redox state in pancreatic β-cells (see Introduction). NADH generated by glycolysis was efficiently reoxidized by highly active mitochondrial shuttles rather than by lactate dehydrogenase in basal conditions in our model (flux JTNAD is considerably higher than JLDH in Figure Figure2A).2A). Steady-state simulation was designed to investigate a role of the redox shuttles in cytoplasmic and mitochondrial events following a step change in the transport rate coefficient for NADH shuttles at several glucose levels (Figure (Figure66).
The result of this simulation showed that [Ca2+]c quickly reached saturation with increased shuttle activity from the basal level. However, a decrease in transport rate coefficient (TNADH) resulted in a decrease in ATP/ADP ratio and in [Ca2+]c from the apparent threshold (Figure (Figure6A).6A). Failure to activate NADH-NAD+ transport between the cytosol and mitochondria does not significantly alter the mitochondrial [NAD+]m/[NADH]m ratio but significantly increases the cytosolic [NADH]c/[NAD+]c ratio, since [NADH]c increases quickly with decreased TNADH (Figure (Figure6B).6B). Increasing the [NADH]c/[NAD+]c ratio led to an acceleration of lactate production from pyruvate (Equation 9, Appendix). The corresponding reduction in pyruvate concentration decreased mitochondrial ATP production, leading to a decreased ATP/ADP ratio and decreased [Ca2+]c. However, glucose consumption did not decrease significantly, because excess glucose influx was shifted to accelerated lactate production in the cytoplasm, as a consequence of the sharp increase in the cytoplasmic NADH concentration (Figure (Figure6).6). This in silico study shows that the cytosolic NADH transport into the mitochondria can be a key regulator of cytosolic NADH/NAD+ and lactate production in pancreatic β-cells only when its influx rate was considerably decreased from basal conditions. Similar results were obtained in a simulation of the shuttle system in cardiomyocytes .
The model simulations on the role of shuttles are also similar to published data. For example, in islets deficient in glycerol-3-phosphate dehydrogenase (GPDH) (which would effectively eliminate the glycerol-phosphate shuttle), when the malate-aspartate shuttle was blocked by inhibiting aspartate aminotransferase by aminooxyacetate, glucose-induced increases in cellular ATP content were impaired and insulin secretion was eliminated, whereas the glycolytic flux remained unchanged . However, studies of this mGPDH-/- mouse model also show that neither absence of the glycerol-phosphate shuttle (in mGPDH-/- islets) nor suppression of the malate-aspartate shuttle alone (in wild-type islets) altered ATP synthesis or GSIS . Our simulations (Figure (Figure6)6) suggest an explanation for these interesting results. Decreased shuttle activity did not lead to a change in [Ca2+]c, if the TNADH initial value was close to a saturated level, and only significant inhibition of shuttle activity, when both shuttles are blocked, resulted in a decrease of [Ca2+]c, whereas the glucose consumption remained unchanged.
The effects of agents known to increase shuttle activities were also examined. Shuttle agonists increased the average [Ca2+]c in mouse islets in the presence of 12 mM glucose . Adenoviral overexpression of the protein Aralarl, a Ca2+ sensitive member of the malate-aspartate shuttle, in both insulin-secreting INS-IE cells and rat pancreatic islets, enhanced glucose-evoked NAD(P)H autofluorescence, Ψm and insulin secretion. Glucose oxidation was enhanced and lactate production was reduced . These experimental results are in reasonably good agreement with the simulated [Ca2+]c, Ψm and [NADH]m increases and decreased lactate output (JLDH) in Figure Figure66 when TNADH was increased from the basal level.
Interestingly, activity and expression of the key enzymes of NADH shuttles were found to be significantly decreased in fetal rat and pig islets compared with adult islets. This can contribute to the inability of fetal β-cells to secrete insulin robustly in response to glucose . Activity of mGPDH, the key enzyme in the glycerol phosphate shuttle, was reduced in islets from patients with type 2 diabetes , and it follows that decreased levels of NADH shuttle activity could also be a possible contributor to β-cell secretory failure . Our simulations confirm the possibility that substantial decreases in activity of NADH shuttles can result in secretory failure of β-cells. In this case our model also predicts an increased lactate production rate in β-cells.
Calcium signaling is associated with mitochondrial uptake of Ca2+ . Recent islet studies have shown steady state resting mitochondrial [Ca2+]m levels are relatively low. Calibrated resting [Ca2+]m in rat islets was approximately 250 nM at 3 mM external glucose. Increasing glucose to 16 mM resulted in a rise of cytoplasmic and mitochondrial [Ca2+] above its resting level [48,50]. The simulated increase of [Ca2+]m versus glucose concentration for steady-state (Figure (Figure2E)2E) was in agreement with these experimental data.
We explored the parameter space with respect to the mitochondrial Na+/Ca2+ antiporter rate to simulate a change of [Ca2+]m. We simulated inhibition and activation of the Na+/Ca2+ exchanger velocity at medium glucose level that resulted in changes in [Ca2+]m at steady-state (Figure (Figure7).7). We analyzed several aspects of the possible role of [Ca2+]c in regulating mechanisms comparing of the experimental evidences and our simulations.
Mitochondrial Ca2+controls the key rate-limiting steps in the TCA cycle through activation of pyruvate dehydrogenase and at least two TCA cycle enzymes: isocitrate dehydrogenase and α-ketoglutarate dehydrogenase (reviewed in ) and F1F0 ATPase . A control hypothesis emerged from this discovery . According to this hypothesis an increase in glucose concentration is accompanied by a rise in cytoplasmic Ca2+, and the subsequent effect of matrix Ca2+ on the TCA cycle increases the supply of reducing equivalents (NADH, FADH2) leading to a "push" of electrons through the respiratory chain. This accelerates ATP production by generating more proton motive force, leading to increased [NADH]m and stimulating oxidative phosphorylation [9,81].
However, a dominant role of [Ca2+]m in a control of β-cell oxidative metabolism under physiological conditions is questionable. The initial mitochondrial NADH and Ψm response precedes increased cytoplasmic and mitochondrial Ca2+ in response to a sharp increase in glucose [3,20,21,50]. In addition, respiration rate and insulin secretion can initially follow the glycolytic rate, which is determined by glucokinase activity, rather than by [Ca2+]m increase (see ). For this reason, it was proposed that Ca2+ is more involved in the maintenance rather than in the initiation of glucose metabolism-secretion coupling. At high glucose, glycolysis may become sufficiently fast such that pyruvate oxidation becomes rate limiting in the formation of ATP or glucose-derived intermediates [41,52]. Under this limitation, positive regulation of mitochondrial metabolism can require additional Ca2+ activation for the synthesis of ATP and other coupling factors (see ).
Because [Ca2+]m can be regulated by mitochondrial Na+/Ca2+ exchanger activity, inhibition of the Na+/Ca2+ exchanger was suggested as a possible target to increase [Ca2+]m and thereby improve insulin secretion in type 2 diabetes. An inhibitor of the exchanger (CGP37157) was shown to prolong mitochondrial Ca2+ signals and increase insulin secretion . (However, this inhibitor could block plasma membrane Ca2+ channels at high concentrations ).
Our model allows an evaluation of the influence of [Ca2+]m changes on GSIS. The results of simulations (Figure (Figure7)7) showed that increasing mitochondrial [Ca2+]m by inhibiting the Na+/Ca2+antiporter did not initially lead to any changes in mitochondrial fluxes or the corresponding increase in the ATP/ADP ratio and [Ca2+]c at all glucose levels evaluated. The initial decrease of [Ca2+]m due to an increased maximal velocity of Na+/Ca2+ also did not lead to significant changes in the ATP/ADP ratio and [Ca2+]c (Figure (Figure7).7). We found that the reason for such insensitivity to [Ca2+]m was that [Ca2+]m was above the threshold for activation of mitochondrial processes even at basal conditions. Our simulations showed that the respiration rate and insulin secretion may follow the glycolytic rate at physiological conditions, rather than an increase in [Ca2+]m (see above).
However, a large decrease in [Ca2+]m due to a large increase in the maximal velocity of Na+/Ca2+ exchange led to an inhibition of ATP production and a decreased ATP/ADP ratio and [Ca2+]c (Figure (Figure7).7). Our simulations also suggest that an increase in GSIS can be expected following increased [Ca2+]m, for example, following inhibition of the Na+/Ca2+ exchanger but only if the initial [Ca2+]m is so low as to limit mitochondrial reactions. The effect of [Ca2+]m on activation of mitochondrial processes in β-cells should therefore be further tested specifically under hyperglycemic and hyperlipidemic conditions.
Physiological influx of Ca2+ into the mitochondrion can cause a measurable concurrent mitochondrial depolarization . Ca2+ cycling in mitochondria reflects the uptake of Ca2+ electrogenically coupled to the efflux of Ca2+ in exchange for protons or sodium ions. This effectively results in uncoupling and would be included in proton leak measurements (Equation 18, Appendix). A high mitochondrial Ca2+ influx and efflux equivalent to an increased proton leak would cause a fall in membrane potential. This would shut down ATP production by the F1F0 ATPase, a process referred to as "short circuiting" in the several mathematical models for β-cell mitochondria in [29,30,85]. In these models the uptake of Ca2+ by β-cell mitochondria suppressed the rate of production of ATP via oxidative phosphorylation.
However, the conclusion that mitochondrial Ca2+ cycling leads to energy dissipation in the pancreatic β-cell is not supported by experimental data which actually favors the opposite idea, that the primary role of mitochondrial Ca2+ is the stimulation of oxidative phosphorylation . The contribution of Ca2+ cycling to proton leak was estimated to be only about 1% of the state 3 rate [9,86]. Mitochondrial Ca2+ cycling also does not lead to marked energy dissipation in the heart . Recordings in mouse islet cells revealed no effect of the inhibitor of the Na+/Ca2+ exchanger (CGP-37157) on the hyperpolarized Ψm under glucose-stimulated conditions  suggesting that mitochondrial Ca2+ cycling likely does not make a major contribution to energy dissipation. Our simulation shows increased ATP production coincidentally with increased [Ca2+]c does not decrease Ψm (Figure (Figure2).2). The contribution of Ca2+ cycling to the fluxes that result in a dissipation of Ψm was less than 1% in our model. This is due to low activity of the Ca2+ uniporter and Na+/Ca2+ exchanger that we have employed to approximate the observed delay between the oscillations [Ca2+]c and [Ca2+]m in β-cells in vivo (see Appendix). The overestimation of the role of Ca2+ fluxes in the dissipation of Ψm in β-cell models [28-30,85] may have arisen from the overestimation of mitochondrial Ca2+ fluxes taken from data obtained on isolated mitochondria.
The role of mitochondrial Ca2+ handling in the regulation of cytoplasmic Ca2+ concentration has been emphasized in several cell types . Mitochondria are an important storage component for Ca2+ handling in cardiomyocytes, where fast and large cytoplasmic Ca2+ changes define cardiac excitation-contraction coupling [87,88]. However, an estimation of the subcellular compartmental volumes of cardiomyocytes gave 58.5% cytosolic and 36% mitochondrial volumes, respectively . On the other hand, β-cell mitochondrial volumes ranged from only 4% to 8% per cell [90,91]. This greatly limits the degree to which mitochondria can regulate cytoplasmic Ca2+ concentration in β-cells.
Ca2+ influx through L-type Ca2+ channels and efflux through plasma membrane pumps, along with endoplasmic reticulum (ER) stores are the principal regulators of β-cell cytoplasmic Ca2+ homeostasis [4,5,10,92,93]. A robust mitochondrial Ca2+ pool was not a necessary component of our previous model examining regulation of cytoplasmic Ca2+ homeostasis (see [5,26]). The mitochondrial Ca2+ pool is significantly smaller compared with the ER free Ca2+ pool, because the ER volume can be up to 20% of total β-cell volume  and its free Ca2+ concentration can reach several hundred mM (see Refs. [4,5,92]). Given the smaller mitochondrial volume, our simulations did not show significant Ca2+ fluxes between cytoplasm and mitochondria (see above). This raises the question as to whether β-cell mitochondrial Ca2+ handling can play a significant role in the regulation of β-cell [Ca2+]c under physiological conditions.
Reduction in mitochondrial metabolism and/or cellular content (number or volume) can in principle underlie progression to the decreased insulin secretion typical of Type 2 diabetes [35,37,94]. For this reason, we employed our model to simulate the effect of changes in mitochondrial functional activity and/or content.
We varied the maximal rate of ETC proton pumping (Vme) (Equation 5, Appendix) to simulate of the results of experiments where the electron transport chain activity was decreased (Figure (Figure8).8). As expected, this simulation showed an increased [NADH]m with decreased Vme, however, Ψm, ATP production, respiration rate, ATP/ADP ratio and [Ca2+]c were also decreased. The glucokinase reaction rate was not changed because the rate of pyruvate reduction by LDH was increased as a consequence of increased [Pyr]. This was due to increased [Pyr] and the part of glycolytic substrate that was not used for ATP synthesis was removed by increased lactate production (not shown).
These simulations correspond to data [95-97] obtained following inhibition of transcription of mitochondrial DNA by ethidium bromide (EtBr), leading to a reduced expression of the mitochondrial electron transport system. After EtBr treatment (in the INS-1 β-cell line), Ψmfailed to hyperpolarize in response to glucose and ATP production and insulin secretion were significantly decreased . Noda et al.  found EtBr caused increased NADH accumulation and lactate production in β HC9 cells, along with decreased [Ca2+]c and ATP/ADP ratio. In contrast, glucose utilization was only insignificantly decreased.
We varied also the maximal F1F0 ATPase activity (Vmph) (Appendix, Equation 7) to simulate of the results of experiments where this activity was changed. Suppression of F1F0 ATPase activity resulted in decreased of ATP synthesis in mitochondria. This decreased [ATP]c and increased [ADP]c leading to a decreased ATP/ADP ratio and decreased [Ca2+]c (Figure (Figure9).9). Decreased proton flux rates through the F1F0 ATPase led to an increase in Ψm. According to Equation 5C an increase in Ψm decreases ETC activity. This decreased [NADH]m consumption lead to decreased respiration rate (measured as O2 consumption) (simulations not shown) and insignificant [NADH]m accumulation. Interestingly, here as well as in first case, the rate of glucokinase reaction was not changed because the NADH not utilized for ATP production was expended on increased leak due to increased Ψm (not shown).
These simulations correspond to the data obtained on BHE/cdb rats, which have a mutation in ATP synthase that limits ATP production and leads to development of mild diabetes [98,99]. BHE/cdb rat islets showed reduced responsiveness to glucose stimulation and ATP content was lower than in control islets . The authors suggested that Ψm is increased in BHE/cdb rat islets due to increased oxygen radical formation . GSIS was reduced, but could not be attributed to changes in glucokinase activity or islet glucose uptake .
To simulate changes in mitochondrial activity or content we varied the mitochondrial volume and corresponding maximum rate of the mitochondrial reaction fluxes in the model (Figure (Figure10).10). The simulation showed that decreased maximal rates of all intramitochondrial processes from the basal level, corresponding to decreased mitochondrial activity or content per cell volume, resulted in increased pyruvate and [NADH]m concentrations, coincidentally with a decreased Ψm, ATP/ADP ratio and [Ca2+]c (Figure 10A,B). The glucose consumption rate decreased moderately since the increase in [PYR] led to an increase in lactate production (Figure 10C) that together with a decreased rate of pyruvate decarboxylation in mitochondria led to an insignificant change in glycolytic flux.
These model simulations for decreased mitochondrial activity or content are in accordance with experimental data. For example, a pancreatic β-cell mouse model for mitochondrial diabetes induced by tissue-specific disruption of the nuclear gene encoding the mitochondrial transcription factor A (Tfam) displayed severe mtDNA depletion, deficient oxidative phosphorylation and abnormally enlarged islet mitochondria . Tfam is essential for mtDNA expression and maintenance, β-cell stimulus-secretion coupling in isolated islets from tfam -/- mice showed reduced hyperpolarization of the mitochondrial membrane potential, impaired Ca2+-signaling and lowered GSIS . Similarly, β-cell specific disruption of Tfam led to 50% reduction in mRNA levels for the mitochondrially encoded nd1 gene, a subunit of the NADH dehydrogenase comprising complex I of the mitochondrial respiratory chain. As a consequence, total cellular ATP concentration was drastically decreased by 75%, and glucose failed to augment cytosolic ATP, explaining the blunted GSIS .
Simulation of increased mitochondrial content above basal levels leads to an initial increase in [Ca2+]c at the tested glucose levels. This initial increase of [Ca2+]c is due to increased ATP production (and an increase in ATP/ADP ratio) as a result of decreased lactate production (see Figure 10C). Interestingly, our simulation shows decreased [Ca2+]c following increased mitochondrial content, due to decreased ATP production (with decreased ATP/ADP ratio) at low and medium glucose levels. This is a result of increased total leak activity that occurs in the model as a consequence of increased mitochondrial content. According to the simulation, lactate production was higher and relative leak was lower (compared with electron transport rate) in basal conditions at high glucose content. For this reason [Ca2+]c did not decrease with increased mitochondrial activity or content at high glucose level ([Glu] = 20 mM) (Figure 10A).
While our simulations support the obvious result that a decrease in mitochondrial function (or content) leads to decreased ATP production, ATP/ADP ratio and [Ca2+]c response to glucose, together giving decreased glucose sensitivity (Figure (Figure8,8, ,9,9, ,10),10), the simulations also show that subtle variation in mitochondrial function or content could underlie β-cell defects in type 2 diabetes (see [10,37,103]. On the other hand, increased mitochondrial content can, in theory, initially increase β-cell sensitivity to glucose (Figure (Figure10)10) especially if their initial content was decreased in comparison with basal levels. This supports the idea that an increase in mitochondrial content can be a possible target for treatment of type 2 diabetes.
Oscillations in cytoplasmic and mitochondrial metabolism, membrane potential, intracellular and mitochondrial Ca2+ due to increased glucose concentrations has been described as a specific characteristic of glucose signaling in the β-cell [3,20-22,104]. The source of these oscillations and the mechanism of orchestration are not clearly understood and may reflect multiple processes [5,23,26,105]. Particularly, [Ca2+]c oscillations can be the driving force for other oscillations in pancreatic β-cells (see [26,106]). For this reason, the aim of this section is testing the hypothesis that independent [Ca2+]c oscillations can create coupled oscillations in metabolic processes.
We previously developed a mathematical model describing β-cell ion regulation that shows how [Ca2+]c oscillations can be independently established from mitochondrial processes when the ATP/ADP ratio achieved some threshold leading to initial depolarization of the plasma membrane . However, for ease of use here we employed a simplified mathematical model that created a periodically varied [Ca2+]c in the cytoplasm (Equations 24-30, Appendix) (Figure 11A). Using this simpler model we simulated the characteristic shape of slow [Ca2+]c oscillations (with a period of one minute and longer), where [Ca2+]c increased sharply and the quiescent time was longer than the period with increased [Ca2+]c (the experimental examples of such oscillations in β-cells are shown [5,23,106]). Our model simulated the corresponding changes in several β-cell mitochondrial and metabolic processes (Figure (Figure11).11). This simulation showed that the metabolic and membrane variables in the cytoplasm and mitochondrial matrix can display oscillatory behavior when [Ca2+]c oscillated independently.
A mechanism involving periodic nucleotide concentrations links [Ca2+]c changes to activation of metabolic oscillations. Each [Ca2+]c increase during oscillations leads to increased cytoplasmic ATP consumption (Equation 20, Appendix). This decreases [ATP]c and increases [ADP]c (Figure 11B) leading to a decreased ATP/ADP ratio. This resulted in an amplification of ATP synthesis by the mitochondrial F1F0 ATPase (Equation 7A). Increasing the rate of proton flux through F1F0 ATPase led to a decrease in Ψm (Figure 11B). According to Equation 5C (Appendix) a decrease in Ψm increases ETC activity. This enhanced [NADH]m consumption and the respiration rate, measured as O2 consumption. For this reason, the electron transport and respiration rates were substantially in phase with [Ca2+]c oscillations (Figure 11C). The oscillations in the glucose consumption rate determined by glucokinase (Equation 1, Appendix), that follows the [ATP]c changes, were out of phase with [Ca2+]c oscillations (Figure 11C).
The experimental data are in accordance with our model simulations. For example, increased [Ca2+]c in MIN6 cells at constant glucose levels caused a fall in the ATP/ADP ratio as inferred from an experiment tracking luciferase-generated photons in transfected cells expressing luciferase . Glucose-induced NAD(P)H and [Ca2+]c slow oscillations were measured simultaneously in mouse pancreatic islets, revealing that NAD(P)H oscillations were small and preceded those of calcium by about 0.1 of a period . In our model the mitochondrial NADH peak concentrations also slightly preceded [Ca2+]c periodic maxima (see Figures 11A,B). The delay (about 9 sec) was about 0.1 of a period of simulated [Ca2+]c oscillations.
Glucose-induced [Ca2+]c and Ψm slow oscillations have been reported [20-22,104], and measured simultaneously in mouse pancreatic islets . The results (Figure (Figure44 from Ref. ) were similar to the dependence that we calculated (Figures 11A,B).
Oscillations in islet oxygen and glucose consumption have also been recorded [95,107]. The glucose consumption rate was out of phase with slow [Ca2+] oscillations, and oxygen consumption rate and [Ca2+]c changes were approximately in phase. The mechanism of this phenomenon is as yet unknown. However, these data are in accord with our model simulation (Figure 11A,C) and could be explained by decreased ATP and increased free ADP concentrations with increased [Ca2+]c during the appropriate phase of slow oscillations.
Mechanisms other than [ADP]c changes can cause variations in [Ca2+]m. Because the changes of Ψm were modest while [Ca2+]c varied with significant amplitude during oscillations, an influx of Ca2+ into mitochondria is determined mainly by a change in [Ca2+]c while Ca2+ efflux depended predominantly on [Ca2+]m changes in our model. This led to some delay in [Ca2+]m oscillations in comparison with [Ca2+]c (14 sec in Figure 11A). The value of this delay was used here to determine maximal rates of Ca2+ fluxes (PCa in Equation 10 and VNC in Equation 11, Appendix).
These results of simulations were consistent with observations demonstrating that oscillations in mitochondrial Ca2+ were in response to glucose elevations, presumably tracking oscillations in [Ca2+]c [50,82,108,109]. Periodic oscillations in [Ca2+]m followed [Ca2+]c oscillations with a delay of approximately 14 sec [50,109]. A similar delay between the maxima of [Ca2+]c and [Ca2+]m (14 sec) was simulated in Figure 11A.
From these results, our simulations at least partially confirm a possible cycle of events previously suggested whereby increased [Ca2+]c during oscillations results in a decrease in ATP/ADP ratio due to increased ATP consumption [13,50]. In the next phase, when [Ca2+]cdecreases, ATP production outweighs ATP consumption leading to an increasing ATP/ADP ratio. However, this pathway is unlikely to be a possible pacemaker mechanism for [Ca2+]c oscillations as has been proposed [13,50]). Rather this could serve as a mechanism where cyclic changes in ATP/ADP ratio are determined by independent cytoplasmic Ca2+ oscillations.
However, they have a typical dynamic that is intrinsic to two successive components where [Ca2+]m follows [Ca2+]c with a particular delay (see above). Our dynamic simulations also clearly show that the independent [Ca2+]c oscillations lead to simulation of [Ca2+]m, Ψm, mitochondrial NADH and respiration oscillations that were similar to experimental observations. These experimental data and our simulations suggest that independent [Ca2+]c oscillations can be a pacemaker in the generation of oscillations of mitochondrial and cytoplasmic parameters in β-cells.
In most cells mitochondria represent the main source of the physiological production of reactive oxygen species (ROS), which may be a byproduct of ETC function [9,110], while recent evidence suggests that NADPH oxidase-dependent generation of ROS in pancreatic β-cells  and ROS generation from sulfhydryl formation in proinsulin biosynthesis  are also potentially important potential sources of ROS generation.
ROS production in mitochondria depends upon the redox state of the ETC complexes, since the ETC carriers in a reduced state have the property of donating electrons to oxygen . The redox state of the ETC complexes and, consequently, the rate of superoxide production also highly depend on Ψm. The increased Ψm (above about 160 mV) decreases electron transport capability leading to a reduced state of the carriers and sharply increases ROS production [110,113]. A pronounced increase in Ψm also augmented ROS production in pancreatic β-cells [36,70].
Interestingly, β-cells have relatively low levels of free radical detoxifying and redox-regulating enzymes such as superoxide dismutase, glutathione peroxidase, catalase and thioredoxin [114,115]. The reasons for this are unclear, although it has been suggested ROS may also subserve a signaling function . However, despite the reported low expression of detoxifying and redox-regulating enzymes in β-cells, antioxidant systems appear to be sufficient to prevent acute oxidative damage under normal physiological conditions [114,115]. We can explain this intriguing property of β-cells based on our simulations results, which show (Figures (Figures22 and and3)3) that β-cells usually work at a relatively lower Ψm(< 160 mV) compared with other types of cells (see above). This leads to decreased ROS production and, consequently, β-cells require a decreased level of detoxifying and redox-regulating enzymes, at least under "normal" conditions, as has been reported.
However, a decreased level of detoxifying and redox-regulating enzymes may also contribute to increased sensitivity of β-cell mitochondria to damage when ROS production increases above normal levels. For example, as was pointed out above, increased respiratory activity, decreased UCP content, decreased LDH activity or high glucose levels could lead to increased Ψm and subsequently to increased ROS production. The imbalance in ROS production versus protection can occur in β-cells more easily than in other types of cells because of decreased level of detoxifying and redox-regulating enzymes.
For these reasons, an overload by metabolic secretagogues such as glucose, causing increased insulin secretion following increase in Ψm and ATP/ADP ratio can also induce increased oxidative stress as a result of elevated ROS production and decreased ability for defense. This could quickly result in ROS-induced reduction in mitochondrial function and cellular content, decreased ATP synthesis, dysregulated calcium homeostasis and ultimately cell death. Decreased mitochondrial functional activity or content could lead to decreased glucose sensitivity and consequently can lead to type 2 diabetes [36,37,117].
The integrated mathematical model of β-cell energetics constructed here has allowed us to reproduce and suggest explanations for experimental relationships among Ψm, respiration, NADH, mitochondrial and cytoplasmic Ca2+ and ATP/ADP ratio and other parameters under various conditions. Our analysis shows that the control structure of β-cell energetics is determined predominantly by the need to regulate the ATP/ADP ratio to keep sensitivity to glucose within the physiological range.
Special features of glucose metabolism in pancreatic β-cells are central to enable this physiological role. Three of these characteristics were emphasized in our model: the glucose phosphorylation by the high-Km glucokinase, which is rate-limiting for glucose-sensitive metabolism and determines the glucose dependency curves of many processes in the β-cell; the remarkably low activity of LDH; and the presence of effective hydrogen shuttles to allow virtually quantitative oxidation of glycolytic NADH (see above). Our analysis is both influenced by and supports these proposals.
We found that the adaptive mechanisms in mitochondria can include decreased respiratory activity compared with F1F0 ATPase activity. This mechanism leads to mitochondrial function in the region with decreased Ψm where variations in Ψm should result in high sensitivity to physiological glucose levels and low ROS production. We found also that proton leak allows fine adjustment and exerts a higher level of control by shifting the ATP/ADP ratio and [Ca2+]c sensitivity to glucose.
This comprehensive model of β-cell energetics permits testing of bioenergetic control hypotheses, provides a basis for further refinement of such steps as oxidative phosphorylation, mitochondrial Ca2+ handling and allow for future incorporation of other biochemical pathways.
The reaction network of the model is shown in Figure Figure1.1. Due to the large number of reactions under consideration, the mathematical description is simplified as much as possible and in some cases intermediate metabolites were lumped together. To help simplify the mitochondrial currents and variables, we represented respiration and ATP synthesis rates in terms of effective proton fluxes.
The parameters for the model are given in Tables Tables22 and and33 along with references to the specific original studies. However, to achieve a closer correspondence with metabolism in intact β-cells or islets we have modified some of these parameters to more closely reflect the complex experimental measurements in vivo, where indicated.
When the concentration of glucose rises in the extracellular medium, the β-cells take up glucose by means of a glucose transporter, GLUT2 and/or GLUT1. The uptake is rapid and believed to effectively equilibrate the extracellular and intracellular concentrations of the sugar . The initial phosphorylation of glucose to glucose 6-phosphate (G6P) is catalyzed by glucokinase (see Introduction).
In the absence of externally added pyruvate or other metabolites, β-cells are glycolytic. For the breakdown of each G6P molecule, the overall reaction yielded two pyruvate molecules, the phosphorylation of three ADP molecules to form ATP and the reduction of two NAD+ molecules to NADH. Glucose phosphorylation is the key limiting step for the steady-state rate of glycolysis in β-cells. However, the overall rate of ATP synthesis at high workloads can also be limited by glycolysis due to a decreased cytoplasmic NAD+ availability for glyceraldehyde 3-phosphate dehydrogenase [18,118].
We simulated the phosphorylation of glucose and glycolytic flux by equations for the limiting enzymes. The glucokinase reaction was modeled as previously :
where Jglu is the rate of the glucokinase reaction, Vmglu is the maximum rate of glucose consumption, [Glu] is extracellular D-glucose concentration, [ATP]c is the cytoplasmic ATP concentration, KmATP is the Michaelis-Menten constant, Kmgl is the glucose concentration at which the reaction is half maximum, and hgl is the Hill coefficient.
We lumped several glycolytic reactions and simulated only the NAD+-dependent step. The flux through glyceraldehydes 3-phosphate dehydrogenase (Jgpd) was assumed to follow irreversible Michaelis-Menten kinetics, where a dependence on [NAD+]c/[NADH]c was written as suggested previously :
where VmGPD is the maximal glyceraldehyde3-phosphate dehydrogenase activity, [G3P] is the cytoplasmic glyceraldehydes 3-phosphate concentration, KmGPD is the Michaelis-Menten constant, [NAD+]c and [NADH]c are the cytoplasmic NAD+ and NADH concentrations, and KGNc is the [NAD+]c/[NADH]c ratio that gives half maximal glycolytic flux and NADH production.
The flux catalyzed by lactate dehydrogenase (JLDH) was described :
where VmLDH is the maximal lactate dehydrogenase activity, [Pyr] is the cytoplasmic pyruvate concentration, KmLD is the Michaelis-Menten constant, KLNC is the [NADH]c/[NAD+]c ratio that gives half lactate and NAD+ production. Basal level of VmLDH was evaluated (Table (Table3)3) on the basis of the proposal that the rate of lactate output is approximately 5% of the rate of glucose consumption at [Glu] = 8 mM and other coefficients were at their basal level.
The cytoplasmic pyruvate concentration is the key factor for mitochondrial ATP synthesis. After being transported into the matrix of mitochondria, pyruvate is decarboxylated in a process catalyzed by the multienzyme complex pyruvate dehydrogenase (PDH), thus forming the intermediate acetyl-coenzyme, or carboxylated by pyruvate carboxylase to form oxaloacetate (anaplerotic pathway) [120,121]. As in a previous model  we assumed that the citric acid dehydrogenase rate is proportional to the reaction rate of the PDH reaction. We assigned PDH a commanding role in the regulation of flux of glycolytic metabolites into the TCA cycle. In this step NAD+ is also reduced to NADH, and PDH activity is regulated by the availability of its coenzyme NAD+, i.e. its activity is decreased when high ratios of NADH/NAD+ prevail . Calcium also activates pyruvate dehydrogenase . We described the rate of pyruvate consumption in mitochondria (JPYR) as the product of several regulated factors:
VmPDH represents the maximal PDH activity, fPYR is the pyruvate kinetic factor, [NAD+]m and [NADH]m are the mitochondrial NAD+ and NADH concentrations, fPNAD is the [NAD+]m/[NADH]m kinetic factors and FPCa is the calcium activity factor, and [Ca +]m is the mitochondrial Ca + concentration.
Equations for fPYR and fPNAD are based on the Equation 18 from , KMPYR is a phenomenological Michaelis-Menten constant, KPNm is the ([NAD+]m/[NADH]m ratio that gives half maximal NADH production. Equation 18 from  was used for FPCa, where u1, u2 and KCam are the parameters for fraction of activated PDH.
Ψm is maintained primarily by the action of respiration-driven proton pumps in the electron transport chain that use the energy contained in NADH and FADH2 to pump hydrogen ions (H+) across the mitochondrial membrane out of the mitochondrial matrix. This process depends on NADH, voltage and oxygen concentration.
As protons are pumped across the inner mitochondrial membrane oxygen is consumed by the ETC. Thus, the respiration-driven proton flux is linked to O2 consumption. Both NADH and FADH2 are electron donors, but in our model NADH is the primary and dominant donor. Therefore we were able to omit a contribution from FADH2. We described respiration in terms of an effective-driven proton flux. We based the equation on the conceptual model  where the steady-state flux of electrons through the ETC was represented as a product of several factors:
Jhres is the rate of proton pumping through ETC, Vme is the rate at optimal conditions. FDe is the donor kinetic factor, fAe is the acceptor kinetic factors and, FTe is the thermodynamic potential factor.
fDe was taken from the model . Note that the apparent affinity of complex 1 for NADH (KmNH in Equation 5A) is very low for in vitro measurements and does not correspond to the measured affinity in vivo . For this reason, we used KmNH = 3 mM (Table (Table3)3) that corresponds to the value found for this reaction in vivo .
Acceptor kinetic factor (fAe) reflects the decrease in electron transport resulting from diminished oxygen availability. In our model oxygen concentration is assumed to be saturated and unlimited.
Thermodynamic potential factor (FTe) determines an inhibition of electron transport with an increase in Ψm. This dependence was evaluated according to  where kAT and kBT are coefficients that were fitted mathematically for mitochondria from muscle cells.
Cellular oxygen consumption rate is an indicator of mitochondrial electron transport. Respiration rate was determined as O2 consumption (JO2). Then the oxygen consumption rate (JO2) can be calculated as
The factor of 0.1 indicates that the full chain from NADH to oxygen is believed to translocate ten protons per oxygen atom consumed where the H+/O ratios for oxidation processes were assumed to be constant .
In an intact cell, under physiological high K+, a contribution of ΔpH to the protonmotive force is usually small [125,126]. For simplicity, we suggest that the F1F0 ATPase uses primarily the mitochondrial membrane potential to generate ATP from ADP and Pi by allowing H+ ions to flow into the mitochondria. Some experimental evidence of this has been reported . The rate of proton flux through mitochondrial F1F0 ATPase (Jph) was also described as the product of the several specific factors:
Vmph is the rate of the proton flux under optimal conditions. AD is the kinetic factor for free cytoplasmic ADP, where [MgADPf]c is the concentration of free cytoplasmic MgADP, KmADP is the activation rate constant, and hph is the Hill coefficient. AT is the thermodynamic potential factor, where Kph is the membrane potential yielding half maximal ATP production, hp is the Hill coefficient. ACa is a sum of kinetic factors describing the effect of calcium concentration, where KPCam is the [Ca +]m binding constant.
Since ADP concentration in mitochondria cannot be readily measured in living cells, we used the dependence of AD on free MgADP in medium that can be calculated using the Hill equation, where the mitochondrial oxidative phosphorylation rate increases with increased free MgADP concentration. The apparent half-saturated concentration is in the range of 20 - 45 μM and the Hill coefficient is about 2 [128,129]. Note that a supply of ADP from cytoplasm is provided by the action of adenine nucleotide translocator (ANT) on the mitochondrial inner membrane. We did not model this translocator because Equation 7A includes the kinetic properties of ANT .
The thermodynamic potential factor (AT) can be represented as a sigmoid dependence of ATP production on mitochondrial membrane potential . However, in contrast to  we used the data obtained by  for mammalian mitochondria to fit the coefficients for this dependence (Figure (Figure12).12). The solid line shows the model fit where Kmp is 131.4 mV and Hill coefficient is 8. AT voltage dependence saturates when Ψm is above 160 mV. This dependence of ATP production on Ψm was similar for mitochondria from muscle and β-cell . Similar dependence of phosphorylation on Ψm was found for rat liver mitochondria , and that is consistent with experimental findings for ox-heart submitochondrial vesicles . ACa is the again the sum of kinetic factors describing the effect of Ca2+ . Note that Jph increases with the electrical gradient (Ψm), the MgADP concentration in cytosol and with the Ca2+concentration in mitochondria.
The mitochondrial membrane clearly leaks protons, decreasing the energy that can be used to drive ATP synthesis. Up to 20% of the basal metabolic rate may be dissipated in this basal leak, always present in mitochondria . Part of the leak is also due to uncoupling proteins (UCPs) that exist in mitochondria to uncouple oxidative phosphorylation, among other possible functions, and the level of their expression varies. To simulate the effect of UCPs we used a special regulated proton leak coefficient. The dependence of the rate of proton leak back into the inner mitochondrial matrix on Ψm was described as follows :
Jhl is the proton leak across the mitochondrial inner membrane, Plb is the basal leak coefficient, and Plr is the regulated leak coefficient, klp is the membrane potential coefficient.
We evaluated the basal leak coefficient in our model by supposing that the basal leak approaches ~20% of the electron transport rate at 160 mV (see Table Table3).3). We set the equation so that the coefficient Plr equals Plb for basal β-cell simulation (Table (Table3),3), leading to a 2-fold increase of Jh1 (i.e up to 40% of the electron transport rate at 160 mV).
Transport of the reducing equivalents derived from glycolysis in β-cells from the cytosol to mitochondria in exchange for NAD+ is primarily through the glycerol phosphate and the malate-aspartate shuttles [2,17]. These shuttles involve oxidation and reduction reactions and include the enzymatic transporters. However, these processes have not been sufficiently described to be included in a model. Therefore, we employed an empirically derived rate expression , where the shuttle activity is expressed as the function of cytoplasmic and mitochondrial NADH/NAD+ ratios:
where TNADH is the transport rate coefficient, KTNc and KTNm are the affinity coefficients.
Ca2+ influx into the mitochondria is mediated by the Ca2+ uniporter that is regulated by the electrochemical gradient. This Ca2+ uniporter is not likely to be saturable by Ca2+ under physiological conditions and exhibits a half-activation constant for cytoplasmic Ca2+ at concentrations greater than 10 μM . The mitochondrial Ca2+ uniporter was described according to  as
where PCa is the permeability of the Ca2+ uniporter, ZCa is the charge of Ca2+, [Ca2+]c is the cytoplasmic Ca2+ concentration, αm, and αi are the mitochondrial and extramitochondrial activity coefficients.
In most cells, including pancreatic β-cells, the main mechanism of Ca2+ extrusion from the mitochondria is the Na+/Ca2+ exchanger [9,82]. The expression for the Na+/Ca2+ exchanger was described according to  as
JNCa is the Na+ - Ca2+ exchange flux, VmNc is the maximal exchanger velocity, [Na+]c and [Na+]m are the cytoplasmic and mitochondrial Na+ concentration, kNaj is the binding constant for sodium and KCaj is the binding constant for Ca2+.
We set Ca2+ uniporter and Na+/Ca2+ exchanger maximum activity from experimental data obtained for living cells. Delay was observed between changes in the cytosolic Ca + concentration and corresponding changes in [Ca2+]m [50,109]. This lag may reflect a delayed response of the [Ca2+]m on [Ca2+]c changes. We have modeled this delay by setting the maximal activities of the Ca2+ uniporter and the Na+/Ca2+ exchanger at levels that replicate experimental slow [Ca2+]m dynamics (see text).
Finally, based on the above relations we determined the dynamic equations for cytoplasmic [G3P], [Pyr], [NADH]c and free [ADP] and the mitochondrial variables: [NADH]m, [Ca +]m and Ψm- All fluxes were expressed as mass fluxes per unit time per unit total cell volume.
The balance equation for the G3P concentration is
where kgpd is the rate constant of G3P consumption in cytoplasm, Vi is the relative cytoplasmic fraction of total cell volume. Coefficient two in the numerator indicates that the breakdown of each glucose molecule yields two G3P molecules.
Pyruvate is the main product of glycolysis. Due to lack of information on pyruvate concentrations at the sub-cellular level in β-cells, we do not differentiate the cytosolic and mitochondrial pyruvate pools in this model. For simplicity, we assume that mitochondrial pyruvate is decarboxylated only in a process catalyzed by PDH. The equation describing [Pyr] change over time is
where Vmmit is the relative mitochondrial matrix fraction of the cell volume.
In this model the TCA cycle is not explicitly modeled. However, it is known that under steady-state conditions, four NADH and one FADH2 are synthesized in the TCA cycle for each pyruvate molecule. To help simplify the mitochondrial variables in our model, we evaluated the reducing equivalents (NADH and FADH2) in terms of H+ flux due to pumping in the ETC. We assumed that 10 and 6 protons are pumped by each NADH and FADH2 oxidation, respectively . This means that one FADH2 can be represented as 0.6 NADH molecule.
We assumed that the NADH production rate in mitochondria is determined by the reaction rate of pyruvate decarboxylation. The mitochondrial NADH concentration is decreased by the action of the ETC, during which NADH is converted to NAD+ and oxygen is consumed. The following equation describes the change in [NADH]m over time
where the stoichiometric coefficient 4.6 in the numerator arises from the number of NADH molecules that are produced from each pyruvate molecule. The coefficient 0.1 indicates that 10 protons are transported by the ETC for each NADH consumed . kNADHm is the rate constant of NADH consumption in mitochondria. The concentration of pyridine nucleotides is assumed to be conserved in the mitochondrial matrix:
where Ntm is the total concentration of pyridine nucleotides in the mitochondrial matrix (Table (Table22).
We assume also that the cytoplasmic concentration of pyridine nucleotides is conserved. The balance equation for the cytoplasmic NADH and NAD+ concentrations is
where Ntc is the total concentrations of pyridine nucleotides in cytoplasm (Table (Table2).2). kNADHc is the rate constant of NADH consumption in the cytoplasm.
Mitochondrial Ψm depends on the respiration module that establishes Ψm, negative inside, and the dissipation processes: the activity of the phosphorylation apparatus, which includes the phosphate carrier, the ATP synthase and the adenine nucleotide translocase and proton-leak reactions. Another portion of Ψm drives Ca2+ out of the matrix via the Na+/Ca2+ exchanger. Ca2+influx into mitochondria mediated by the Ca2+ uniporter is also electrogenic. The final equation for Ψm is
where Cmi is the mitochondrial membrane capacitance. The stoichiometric coefficient 2 in the numerator arises from the number of the charges that flow into mitochondria with one Ca2+. We assumed that Ψm > 0.
The balance equations for [Ca2+]m is described by
where fm is the constant describing the fraction of free Ca2+ in mitochondria.
Several assumptions were made to simplify the ATP regulation model: (1) Due to the rapid export of mitochondrial ATP to the cytosol via the ATP/ADP transporter no limitation for this process was suggested. (2) ADP binding is in equilibrium (i.e. the binding reactions are fast compared to the other processes included in the model of ATP dynamic).
Utilization of ATP in the β-cells is mostly for ion transport, biosynthesis, and secretion. The rate of ATP utilization is a complex function of the concentrations of ATP, [Ca2+]c, glucose and numerous other factors, but almost nothing is known about the quantitative relation between them. Clearly, there is some basal level of ATP consumption in a cell at low glucose and [Ca2+]c concentrations. There is evidence that ATP hydrolyzing processes are accelerated during a glucose-induced increase in [Ca2+]c in β-cells, for example, ATP consumption can be increased due to increased Ca2+ pump activity in the plasma membrane and endoplasmic reticulum [7,10,12,13]. For this reason, two terms for ATP consumption were introduced: basal and Ca2+-dependent (see ). Then, on the basis of the Equation 7 for oxidative phosphorylation we can write the balance equation for [ATP]c and [ADP]c
The factor of 2 in the numerator of Equation 20 (at Jglu) indicates that two ATP molecules are produced during glycolysis from each glucose molecule. The coefficient 0.231 (at Jpf) arises from the ATP/H+ coupling stoichiometry of 3/13 in mammalian mitochondria. This suggested as well the additional transport of one proton from the cytoplasm to the matrix that is associated with the movement of phosphate through the mitochondrial membrane ; KATP is the permanent rate constant of basal ATP consumption, and KATPca is the rate constant of ATP consumption that accelerates as Ca2+ increases.
The general concentration of intracellular nucleotides (A0) was assumed to be constant. The constants KATP and KATP,Ca (Equation 20) were fitted to simulate both the low and high glucose modeled rate of ATP production and ATP/ADP ratio in β-cells (see [5,26]) (Table (Table3).3). The coefficient 0.055 in Equation 22 was calculated as in our previous model [5,26].
The mechanisms of [Ca2+]c and [Na+]c regulation and their interrelationships with other metabolic and ionic fluxes are incompletely understood. For this reason in our model we used their empirical determination as the external parameters, even though we previously developed a mathematical model for [Ca2+]c and [Na ]c regulation in β-cells [5,26]. [Na+]c was proposed as constant in this model (Table (Table22).
The relationship between [Ca2+]c and glucose at steady-state was calculated using the change of ATP/ADP ratio as an intermediate agent by an empirical Hill-type equation.
where [Ca2+]R is the constant [Ca2+] for resting low glucose, kACa is the coefficients, KAD is the [ATP]c/[ADP]c ratio that gives half maximal [Ca2+]c, hCa is the Hill coefficient. The parameters of Equation 23 were set to a value (Table (Table3)3) that produced a reasonable amplification of [Ca2+]c with increased glucose (see Results and Discussion).
We used a simple mathematical model that creates a periodically varied independent [Ca2+]c change in cytoplasm to simulate the oscillations of the mitochondrial parameters in β-cells. The model described in this section is based on a simple model  and used only for simulation of independent Ca2+ oscillations in the cytoplasm of a mean individual cell.
gmVCa is the maximum conductance for IVCa (gmVCa = 20 μS), VP is the plasma membrane potential.
where gmKCa is the maximum conductance for IKCa, KKCa is the half-maximum Ca2+ binding constant for IKCa (gmKCa = 25 μS, KKCa = 0.25 μM).
The differential equation describing time-dependent changes in the plasma membrane potential (Vp) is the current balance equation:
where Cm is the cell membrane capacitance (Cmp = 6158 fF )
The equations for [Ca2+]c dynamics can be written as:
where fi is the fraction of free Ca2+ in cytoplasm, F is Faraday's constant, Vci is the volume of the cytosolic compartment in single cell (0.764 pL ), and ksg is the coefficient of the sequestration rate of [Ca2+]c (ksg = 0.00002 ms-1 ).
For computational purposes we considered the β-cell as an assemblage of mitochondria with similar properties. The units used in the model are time in millisecond (ms), voltage in millivolts (mV), concentration in micromoles/liter (μ, M), flux in μmol ms-1.
A factor (0.31) was used to convert picomoles per islet from metabolism experiments to the cytoplasmic millimolar terms of a single β-cell (as calculated in Ref. ). The mitochondrial protein density for total mitochondria volume is estimated at ~0.3 mg protein/μl, and free water volume in mitochondrial matrix space can be estimated to be 0.24 of total mitochondria volume . In line with these data we used the factor 1.25 to convert the measured nanomoles per milligram mitochondrial protein (nmol mg-1 protein) to mitochondrial matrix space in millimollar terms. In our model the fluxes were specified for the volume of a cell. Multiplying the mitochondrial matrix volume (Table (Table2)2) by its protein density (1.25 g protein ml-1) by unit flux [in nmol (mg protein min-1)] gives the total flux in the mitochondria matrix per unit cell volume.
The general concentration of intramitochondrial and cytoplasmic adenine and purine nucleotides were kept constant during simulations (Table (Table2).2). Model parameters were found by several methods. Specifically, they were obtained from the scientific literature when possible and were also found by fitting specific model equations to experimental data. The third method was to estimate the parameter so that model variable values and time courses closely matched experimental data. Several enzymatic activity values were treated as adjustable parameters, which were adjusted using the reaction stoichiometries to reflect the rate of glucose phosphorylation by glucokinase (Table (Table3).3). Parameter values from Tables Tables22 and and33 were used unless otherwise mentioned. Nine ordinary differential equations (Equations 12-14, 16, 18-20, 29, 30) describe the behavior of [G3P], [Pyr], [NADH]m, [NADH]c, Ψm, [Ca2+]m, [ATP]c, Vp and [Ca2+]c. Coefficients are shown in Tables Tables22 and and3.3. Calculations on Figure Figure1212 and a generation of all figures were performed using "Igor" (IGOR Pro, WaveMetrics, Inc, Lake Oswedo, OR, USA) and Microsoft Excel X.
Simulations were performed as noted previously for an idealized mean individual cell using the software environment from "Virtual Cell" (Fridlyand et al. [5,26]). To calculate the steady-state cellular parameters, the model was allowed to run for at least 10s with no external stimulation. Calculations obtained with the coefficients from Tables Tables22 and and33 have been mentioned in the text as a simulation at basal levels.
This model is available for direct simulation on the website "Virtual Cell" http://www.nrcam.uchc.edu in "MathModel Database" on the "math workspace" in the library "Fridlyand" with the name "GlucoseSensitivity-1" for the general model and with name "GlucoseSensitivity-2" for the general model that also includes independent [Ca2+]c oscillations.
ATP/ADP: ratio of ATP to ADP; ETC: electron transport chain; EtBr: ethidium bromide; GK: glucokinase; Glu: glucose; GSIS: glucose-stimulated insulin secretion; G3P: glyceraldehyde 3-phosphate; KATP: ATP-sensitive K+ channels; LDH: lactate dehydrogenase; PDH: pyruvate dehydrogenase; PYR: pyruvate; ROS: reactive oxygen species; TCA: tricarboxylic acid; Tfam: mitochondrial transcription factor A; UCP2: uncoupling protein 2; Ψm: mitochondrial membrane potential; Subscript c: cytoplasmic compartment; m: mitochondrial compartment.
The authors declare that they have no competing interests.
All authors (LF, LP) contributed equally to the formulation of the model, the estimation of parameters, the biological interpretations and conclusions, and the writing and editing of the manuscript. LF performed construction and simulation of the initial mathematical model. Both authors read and approved the final manuscript.
We thank Dr. N. Tamarina for helpful discussions. This work has been partially supported by National Institute of Diabetes and Digestive and Kidney Diseases Grants DRTC P60DK020595, DK-48494, DK-063493 and a Research Grant from the American Diabetes Association to LP.