Advertisement

A computational model of cardiomyocyte metabolism predicts unique reperfusion protocols capable of reducing cell damage during ischemia/reperfusion

  • Matthias Grass
    Affiliations
    Department of Mechanical Engineering, ETH Zurich, Zurich, Switzerland

    Department of Pathology, Brigham and Women’s Hospital, Boston, Massachusetts, USA

    Program in Human Biology and Translational Medicine, Harvard Medical School, Boston, Massachusetts, USA
    Search for articles by this author
  • Anthony D. McDougal
    Affiliations
    Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
    Search for articles by this author
  • Adriana Blazeski
    Affiliations
    Department of Pathology, Brigham and Women’s Hospital, Boston, Massachusetts, USA

    Program in Human Biology and Translational Medicine, Harvard Medical School, Boston, Massachusetts, USA

    Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
    Search for articles by this author
  • Roger D. Kamm
    Affiliations
    Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA

    Department of Biological Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
    Search for articles by this author
  • Guillermo García-Cardeña
    Correspondence
    For correspondence: Guillermo García-Cardeña; C. Forbes Dewey Jr
    Affiliations
    Department of Pathology, Brigham and Women’s Hospital, Boston, Massachusetts, USA

    Program in Human Biology and Translational Medicine, Harvard Medical School, Boston, Massachusetts, USA

    Cardiovascular Disease Initiative, The Broad Institute of MIT and Harvard, Cambridge, Massachusetts, USA
    Search for articles by this author
  • C. Forbes Dewey Jr.
    Correspondence
    For correspondence: Guillermo García-Cardeña; C. Forbes Dewey Jr
    Affiliations
    Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA

    Department of Biological Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
    Search for articles by this author
Open AccessPublished:February 10, 2022DOI:https://doi.org/10.1016/j.jbc.2022.101693
      If a coronary blood vessel is occluded and the neighboring cardiomyocytes deprived of oxygen, subsequent reperfusion of the ischemic tissue can lead to oxidative damage due to excessive generation of reactive oxygen species. Cardiomyocytes and their mitochondria are the main energy producers and consumers of the heart, and their metabolic changes during ischemia seem to be a key driver of reperfusion injury. Here, we hypothesized that tracking changes in cardiomyocyte metabolism, such as oxygen and ATP concentrations, would help in identifying points of metabolic failure during ischemia and reperfusion. To track some of these changes continuously from the onset of ischemia through reperfusion, we developed a system of differential equations representing the chemical reactions involved in the production and consumption of 67 molecular species. This model was validated and used to identify conditions present during periods of critical transition in ischemia and reperfusion that could lead to oxidative damage. These simulations identified a range of oxygen concentrations that lead to reverse mitochondrial electron transport at complex I of the respiratory chain and a spike in mitochondrial membrane potential, which are key suspects in the generation of reactive oxygen species at the onset of reperfusion. Our model predicts that a short initial reperfusion treatment with reduced oxygen content (5% of physiological levels) could reduce the cellular damage from both of these mechanisms. This model should serve as an open-source platform to test ideas for treatment of the ischemia reperfusion process by following the temporal evolution of molecular concentrations in the cardiomyocyte.

      Keywords

      Abbreviations:

      ODE (ordinary differential equation), RET (reverse electron transport), ROS (reactive oxygen species), TCA (tricarboxylic acid)
      A major cause of death globally, ischemic heart disease, is caused by a narrowing of the coronary arteries over time (
      GBD Causes of Death Collaborators
      Global, regional, and national age-sex-specific mortality for 282 causes of death in 195 countries and territories, 1980-2017: A systematic analysis for the Global Burden of Disease Study 2017.
      ). Tissues become ischemic when blood vessels are not able to provide them with sufficient nutrients and oxygen (
      • Nowbar A.N.
      • Gitto M.
      • Howard J.P.
      • Francis D.P.
      • Al-Lamee R.
      Mortality from ischemic heart disease.
      ). A myocardial infarction caused by a major cessation of coronary blood flow due to arterial dysfunction is a dramatic example of the onset of ischemia. The drop in intracellular oxygen levels in the cardiac muscle cells leads to impaired energy (ATP) production and an inability to sustain contractile function (
      • Ashrafian H.
      • Frenneaux M.P.
      • Opie L.H.
      Metabolic mechanisms in heart failure.
      ). Depending upon the degree and duration of oxygen loss, the affected heart muscle cells can die. The most obvious countermeasure to ischemia is the timely reperfusion of the affected vessel and restoration of normal oxygen levels. The availability of oxygen induces a reactivation of the electron transport chain and production of ATP within the mitochondria.
      Paradoxically, conditions present during ischemia prime the tissue for significant damage upon reperfusion depending on the duration and severity of the ischemic period (
      • Eltzschig H.K.
      • Eckle T.
      Ischemia and reperfusion–from mechanism to translation.
      ,
      • Henry T.D.
      • Archer S.L.
      • Nelson D.
      • Weir E.K.
      • From A.H.
      Postischemic oxygen radical production varies with duration of ischemia.
      ,
      • Kalogeris T.
      • Baines C.P.
      • Krenz M.
      • Korthuis R.J.
      Ischemia/reperfusion.
      ). Excessive amounts of reactive oxygen species (ROS) are produced at various sites in the mitochondria (
      • Hausenloy D.J.
      • Yellon D.M.
      Myocardial ischemia-reperfusion injury: A neglected therapeutic target.
      ,
      • Murphy M.P.
      How mitochondria produce reactive oxygen species.
      ). In the mitochondria's outer membrane, monoamine oxidase (MAO) and cytochrome b5 reductase (Cb5R) can produce ROS (
      • Angelova P.R.
      • Abramov A.Y.
      Functional role of mitochondrial reactive oxygen species in physiology.
      ). In the inner mitochondrial membrane, glycerol-3-phosphate dehydrogenase, linked to the CoQ pool, can contribute to ROS production (
      • Murphy M.P.
      How mitochondria produce reactive oxygen species.
      ,
      • Tretter L.
      • Takacs K.
      • Hegedus V.
      • Adam-Vizi V.
      Characteristics of alpha-glycerophosphate-evoked H2O2 generation in brain mitochondria.
      ). Other enzymes that interact with the mitochondrial NADH pool, such as αKDGH (α-ketoglutarate dehydrogenase) and other enzymes of the tricarboxylic acid (TCA) cycle (
      • Angelova P.R.
      • Abramov A.Y.
      Functional role of mitochondrial reactive oxygen species in physiology.
      ), can contribute to superoxide production in the mitochondrial matrix (
      • Murphy M.P.
      How mitochondria produce reactive oxygen species.
      ,
      • Starkov A.
      An update on the role of mitochondrial α-ketoglutarate dehydrogenase in oxidative stress.
      ). Finally, mitochondrial enzymes not associated with the NADH or CoQ pools, such as cytochrome P450 enzymes, can also generate ROS in the mitochondria (
      • Hanukoglu I.
      Antioxidant protective mechanisms against reactive oxygen species (ROS) generated by mitochondrial P450 systems in steroidogenic cells.
      ), as can NOX4, an NADPH oxidase that produces hydrogen peroxide, although its presence in the mitochondria is debated (
      • Cadenas S.
      ROS and redox signaling in myocardial ischemia-reperfusion injury and cardioprotection.
      ,
      • Hirschhäuser C.
      • Bornbaum J.
      • Reis A.
      • Böhme S.
      • Kaludercic N.
      • Menabó R.
      • Di Lisa F.
      • Boengler K.
      • Shah A.M.
      • Schulz R.
      • Schmidt H.H.H.W.
      NOX4 in mitochondria: Yeast two-hybrid-based interaction with complex I without relevance for basal reactive oxygen species?.
      ,
      • Beretta M.
      • Santos C.X.
      • Molenaar C.
      • Hafstad A.D.
      • Miller C.C.
      • Revazian A.
      • Betteridge K.
      • Schröder K.
      • Streckfuß-Bömeke K.
      • Doroshow J.H.
      • Fleck R.A.
      • Su T.-S.
      • Belousov V.V.
      • Parsons M.
      • Shah A.M.
      Nox4 regulates InsP3 receptor-dependent Ca2+ release into mitochondria to promote cell survival.
      ). Importantly, mitochondrial protein complexes I and III have been identified as the main ROS producers during reperfusion. Chouchani et al. (
      • Chouchani E.T.
      • Pell V.R.
      • Gaude E.
      • Aksentijević D.
      • Sundier S.Y.
      • Robb E.L.
      • Logan A.
      • Nadtochiy S.M.
      • Ord E.N.J.
      • Smith A.C.
      • Eyassu F.
      • Shirley R.
      • Hu C.H.
      • Dare A.J.
      • James A.M.
      • et al.
      Ischaemic accumulation of succinate controls reperfusion injury through mitochondrial ROS.
      ) have shown how the accumulation of the Krebs cycle intermediary succinate during ischemia leads to a reversal of the electron transport chain at complex I and consequently to an excessive superoxide production and oxidative damage. On the other hand, Korge et al. (
      • Korge P.
      • Calmettes G.
      • John S.A.
      • Weiss J.N.
      Reactive oxygen species production induced by pore opening in cardiac mitochondria: The role of complex III.
      ) suggest that reverse electron transport (RET) is only partially responsible for the oxidative damage during reperfusion and highlight the superoxide production at complex III as a key driver of ROS production and cell damage. Other data have been put forward to justify the failure of calcium regulation to be a driving force for dysfunction of the mitochondria (
      • Malyala S.
      • Zhang Y.
      • Strubbe J.O.
      • Bazil J.N.
      Calcium phosphate precipitation inhibits mitochondrial energy metabolism.
      ).
      Although many studies have identified the important role of the mitochondria in ischemia/reperfusion injury (
      • Jassem W.
      • Fuggle S.V.
      • Rela M.
      • Koo D.D.
      • Heaton N.D.
      The role of mitochondria in ischemia/reperfusion injury.
      ,
      • Honda H.M.
      • Korge P.
      • Weiss J.N.
      Mitochondria and ischemia/reperfusion injury.
      ,
      • Lesnefsky E.J.
      • Moghaddas S.
      • Tandler B.
      • Kerner J.
      • Hoppel C.L.
      Mitochondrial dysfunction in cardiac disease: Ischemia–reperfusion, aging, and heart failure.
      ), few have quantitatively shown the mechanisms and conditions present during the initial moments of reperfusion. These initial conditions require a careful assessment of the cumulative effects of ischemia and the time-dependent evolution of molecular components of the cardiomyocyte and its mitochondria. Analytical models should include the time course of restarting blood flow as this may be crucial to the fate of the cardiomyocytes. McDougal and Dewey (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ) created a mathematical model of the anaerobic cardiomyocyte metabolism during the transition from a normal physiological state to one dominated by ischemia. This study extends that model using a biochemical description of the mitochondria that includes some aspects of oxidative phosphorylation and compartmental transport processes, inspired by the work of Beard (
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ). The resulting model predicts the metabolic changes during the critical transition periods between ischemia and reperfusion and, additionally, suggests novel reperfusion protocols to minimize injury from oxidative damage.

      Results

      Static experiments (physiological steady state and sensitivity)

      In order to validate the predictions of our metabolic model, the chemical concentrations of intracellular species related to the glycolytic pathway were computed until a chemical steady state was reached. These values were then compared with experimental results from Kashiwaya et al. (
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      ), who measured glycolytic intermediates in perfused rat hearts (see Fig. 1).
      Figure thumbnail gr1
      Figure 1Validation of glycolytic intermediates. Comparison of normal physiology steady-state concentrations between model predictions (orange) and values obtained experimentally in perfused rat hearts (black) (
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      ).
      Simultaneously, a sensitivity analysis was conducted, where the dependence of chemical concentrations with respect to enzyme activities was computed and quantitatively compared with experimental observations. More specifically, the relative sensitivity of species i with respect to the activity of enzyme j was calculated as the relative change in concentration of species i due to a 1% change in activity of enzyme j. These values were then summed over all species to quantify the influence of a given enzyme/reaction. The summed sensitivity values associated with enzyme j are represented by Sj. Table 1 shows the qualitative distribution of control across the metabolic network and highlights the regulation of ATP levels through multiple mechanisms. The aggregated sensitivity values highlight the importance of the reaction catalyzed by hexokinase for the downstream part of the glycolytic pathway during physiological conditions. With reduced activity of hexokinase, the amount of G6P will be reduced and the cardiac glucose metabolism is expected to slow down. However, Table 1 also shows the distribution of metabolic control among multiple pathways other than glycolysis. While, during ischemia, the sensitivity values are led by the activity of glucose-6-phosphate dehydrogenase, after reperfusion the sensitivity values are very similar to the preischemic values due to the absence of any irreversible processes implemented in our model. The comparison of the sensitivity values reported in Table 1 with experimental studies regarding the distribution of metabolic control in cardiomyocytes such as reported by Kashiwaya et al. (
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      ) further serves as a qualitative validation of our model. For a mathematical explanation of these sensitivity coefficients, please refer to section Simulations/Sensitivity analysis in Computational Methods.
      Table 1Enzymes/reactions exerting the largest overall control
      Perfusion stateEnzyme/reactionPathwaySRank
      PhysiolHKGlycolysis16.541
      ANTATP/ADP Transporter10.732
      ATPaseATP Consumption10.193
      PFKGlycolysis6.934
      G6PDHPentose Phosphate4.435
      IschemiaG6PDHPentose Phosphate17.161
      HKGlycolysis15.582
      PFKGlycolysis5.833
      TPIGlycolysis4.574
      AKEnergy Buffer4.315
      ReperfusionHKGlycolysis16.361
      ANTATP/ADP Transporter10.802
      ATPaseATP Consumption8.293
      PFKGlycolysis7.364
      G6PDHPentose Phosphate4.295
      Values are L1-norm of the relative sensitivities across metabolites.

      Dynamic validations (ischemia and reperfusion transitions)

      To compare model predictions with published experimental results, simulations were carried out using different values for ischemic oxygen concentrations and time durations of the ischemic period. First, results were compared with ischemic canine hearts (
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      ), where creatine phosphate was measured during 36 s of ischemia (left arrow indicates start of ischemia, right arrow indicates end of ischemia) for 100 s in total (see Fig. 2).
      Figure thumbnail gr2
      Figure 2Validation of ischemia and reperfusion transitions in canine hearts. Comparison of cytosolic creatine phosphate concentrations between model predictions (orange) and values obtained experimentally in canine hearts subject to 36 s of ischemia and 64 s of reperfusion. For ischemia simulations, extracellular oxygen concentration was set to 1% of the physiological level.
      Additionally, results were compared with experimental measurements by Clarke et al. (
      • Clarke K.
      • O'Connor A.J.
      • Willis R.J.
      Temporal relation between energy metabolism and myocardial function during ischemia and reperfusion.
      ), where ischemia was introduced in guinea pig hearts for 150 s, and phosphate-related species were measured for the 150 s of ischemia and 150 s of reperfusion (see Fig. 3).
      Figure thumbnail gr3
      Figure 3Validation of ischemia and reperfusion transitions in guinea pig hearts. Comparison of cytosolic phosphate metabolite concentrations between model predictions (orange) and values obtained experimentally in guinea pig hearts subject to 150 s of ischemia (1% of physiological oxygen) and 150 s of reperfusion. The y-axes are normalized to represent % of the total change over 150 s, as reported by Clarke et al. (
      • Clarke K.
      • O'Connor A.J.
      • Willis R.J.
      Temporal relation between energy metabolism and myocardial function during ischemia and reperfusion.
      ) (e.g., Pirel(t)=|Pi(t)Pimin||PimaxPimin|100). ATP:ADP indicates the ratio of ATP to ADP, a common readout to determine the metabolic state of the cardiomyocyte (
      • Maldonado E.N.
      • Lemasters J.J.
      ATP/ADP ratio, the missed connection between mitochondria and the Warburg effect.
      ). A, relative change in cytosolic inorganic phosphate (Pi) in simulations correlates with experimental measurements during ischemia. B, reperfusion simulations of inorganic phosphate (Pi) also show correlation with experiments. Experimental measurements reach peak slightly faster than simulations. C, relative change in cytosolic creatine phosphate (CrP) in simulations correlates with experimental measurements during ischemia. D, similar to phosphate measurements, simulated changes in creatine phosphate (CrP) during reperfusion show correlation with experiments but xperimental measurements reach peak concentrations faster (~50 s) than simulations (~60 s). E, relative change in ratio of cytosolic ATP:ADP in simulations correlates again with experimental measurements. Ratio saturates faster in simulations than in experimental measurements. F, simulated ATP:ADP ratio during reperfusion shows different behavior than experiments. While simulation shows an S-curve with an inflection point at around ~60 s, experiments show a fast saturation of ATP:ADP ratio before stabilizing at ~50%.

      Oxygen dependence of ischemia/reperfusion

      Even if a blood vessel is fully occluded, the neighboring tissue will always have a small residual perfusion of other nonoccluded blood vessels. Thus, the extracellular oxygen concentration will likely not reach absolute 0. Model simulations were carried out at decreased oxygen concentrations between 10% and 0.5% of normal physiological levels, uncovering trends for the cellular mechanisms that correlate with the severity of ischemia (see Fig. 4). The total simulation time was 3000 s, while introducing varying degrees of ischemia at t = 1000 s and reintroducing preischemic oxygen levels at t = 2000 s. Surprisingly, simulation results indicate that the ATP production in cardiomyocytes can meet demand at oxygen levels as low as 10% of preischemic levels for a sustained duration (see Fig. 4A).
      Figure thumbnail gr4
      Figure 4Ischemia/reperfusion simulations for different levels of residual ischemic oxygen concentrations. Oxygen dependence of metabolic measurements during ischemia (first vertical lines; left arrows) and reperfusion (second vertical lines; right arrows). Oxygen levels in legend refer to values in percent of normal physiological level. A, molecular concentration of ATP, (B) reaction rate at complex I of the electron transport chain (i.e., the amounts of NADH and Q are converted to NAD and QH2 per second; see Equation ), (C) mitochondrial membrane potential, (D) reaction rate at complex III (i.e., the amounts of QH2 and cytochrome C (oxidized) are converted to Q and cytochrome C (reduced) per second; see Equation ).
      For ischemic oxygen concentrations lower than 5% of physiological levels, model simulations predict that cellular ATP levels are depleted after approximately 400 s, similar to results reported by McDougal and Dewey (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ) and Ch’en et al. (
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ). Below 5%, energy buffering enzymes such as creatine kinase (CK) and adenylate kinase (AK) initially keep ATP levels at normal physiological levels by using intracellular stores of creatine phosphate, glycogen, and ultimately ADP. Interestingly, the qualitative behavior of the recovery of ATP, i.e., the temporal trajectory, seems to be the same, regardless of ischemic oxygen concentrations, while only the starting point of the ATP production is postponed with decreasing oxygen concentrations (see Fig. 4A).
      After reintroducing preischemic oxygen levels, the ATP stores are already fully recovered after approximately 200 s of reperfusion, regardless of the severity of ischemia. Additionally, in Figure 4B we show how mitochondrial ATP production via oxidative phosphorylation quickly stops at the onset of ischemia due to a reduced activity of protein complexes and consequently a reduced proton-motive force. The anaerobic metabolism takes over the ATP production during the first minutes of ischemia. At approximately 350 s of ischemia ATP levels approach 0 for ischemic oxygen concentrations below 1%. Due to their importance in the scope of ROS production and subsequent oxidative damage, we also tracked the reaction rates at protein complex I and III of the electron transport chain. While they are tightly coupled, the results indicate that the reaction rate at complex III approaches 0 within the first seconds of ischemia, whereas complex I seems to be working at a greatly reduced activity throughout the period of ischemia (see Fig. 4, B and D respectively).
      The model predicts that at ischemic oxygen concentrations above or equal to 5% of the normal physiological level, the overall reaction rate at complex I decreases instantly at the onset of ischemia and also recovers quickly at the onset of reperfusion. However, when reperfusion follows an ischemic period at even lower ischemic oxygen concentrations, the reaction rate at complex I becomes slightly negative before recovering to preischemic levels, indicating that complex I runs in reverse for a short period of time. Again, the simulations show highly nonlinear responses to different ischemic oxygen concentrations at the onset of reperfusion. Additionally, the reaction rate at complex III exhibits a large spike during the first seconds of reperfusion to levels approximately threefold of the preischemic values.
      In Figure 4C, we plot the polarization of the mitochondrial membrane potential during ischemia and reperfusion. Being directly affected by the reactions at complex I and III, the membrane potential also decreases rapidly at the onset of ischemia before asymptotically approaching 55 mV. Upon reperfusion, the membrane potential spikes to a level twice the preischemic value for approximately 2 min before settling at preischemic levels. If ischemic oxygen concentrations stay at 5% or above, the membrane potential exhibits a less drastic decrease at the onset of ischemia, and upon reperfusion the potential does not overshoot. Thus, there appears to be a threshold of oxygen during ischemia that leads to hyperpolarization of the mitochondrial membrane upon reperfusion. Interestingly, the magnitude of the potential spike seems to be independent of the ischemic oxygen levels past a certain threshold, whereas the duration of the elevated membrane potential is inversely related to ischemic oxygen levels.

      A two-step reperfusion protocol

      In order to reduce the sources of reperfusion injury, a two-step reperfusion protocol was simulated (see Figs. 5 and 6). First, ischemia was introduced for 1000 s. Instead of going back to physiological oxygen levels immediately after reperfusion, an intermediate step with oxygen levels not fully at 100% was added for another 1000 s. In Figure 5, the mentioned two-step reperfusion protocol is simulated for different ischemic oxygen concentrations and a constant first-step reperfusion oxygen concentration of 5%. Figure 5B shows that the dominant part of the increase in mitochondrial membrane potential is regained with even as small an oxygen level of reperfusion as 5% of physiological oxygen (i.e., the membrane potential recovers already during the first reperfusion step, regardless of the ischemic oxygen values). The complex I negative reaction rate (reverse electron transport) shows up distinctly at the initiation of the first step of reperfusion, and its duration increases with decreasing ischemic oxygen concentrations (see Fig. 5C). Similar to Figures 5 and 6 shows the influence of the oxygen concentration during the first wave of reperfusion. Here, the ischemic oxygen concentration was set to 0.5% of physiological values and followed by various two-step reperfusion scenarios, with the first reperfusion step ranging from 1% to 100% of physiological O2 levels. The spike in mitochondrial membrane potential was suppressed for cases when 5% to 10% oxygen was utilized for the first reperfusion step (see Fig. 6B). Similarly, Figure 6C shows that reverse electron transport at complex I is minimized for the same range of oxygen levels; and lastly, Figure 6D indicates that this range of oxygen levels results in minimal spiking of the reaction rate at complex III. All of these observations suggest that there is an optimal range for the oxygen concentration during the first step of the reperfusion that minimizes ROS production and consequently reperfusion injury.
      Figure thumbnail gr5
      Figure 5Ischemia/reperfusion simulations for different levels of residual ischemic oxygen concentrations and a two-step reperfusion protocol. Oxygen levels during ischemia ranged from 0.2% to 10% of physiological levels. At 2000 s, oxygen levels were adjusted to 5% of physiological levels (first step of reperfusion) and at 3000 s, oxygen levels were returned to physiological levels (second step of reperfusion). Traces of (A) oxygen concentrations, (B) mitochondrial membrane potential, (C) reaction rate at complex I, and (D) reaction rate at complex III are color-coded based on the initial level of ischemia.
      Figure thumbnail gr6
      Figure 6Ischemia/reperfusion simulations for different levels of residual oxygen concentrations during the first step of a two-step reperfusion protocol. In each simulation, the oxygen level during ischemia was maintained at 0.5% of physiological levels. At 2000 s, the oxygen level was adjusted to a concentration ranging from 1 to 100% of physiological levels (first step of reperfusion) and at 3000 s, oxygen levels were returned to physiological levels (second step of reperfusion). Traces of (A) oxygen concentrations, (B) mitochondrial membrane potential, (C) reaction rate at complex I, and (D) reaction rate at complex III are color-coded based on the oxygen level during the first step of reperfusion.

      Discussion

      Model validation

      The model developed herein has been validated with experimental values obtained during normal physiological and pathological conditions and from a variety of species. Figure 1 shows that the extensions introduced to the McDougal model (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ) result in similar steady-state values for cytosolic metabolites as measured in perfused rat hearts by Kashiwaya et al. (
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      ). Additionally, model predictions for ischemia and reperfusion conditions correlate well with experimental values measured by Wu et al. (
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      ) and Clarke et al. (
      • Clarke K.
      • O'Connor A.J.
      • Willis R.J.
      Temporal relation between energy metabolism and myocardial function during ischemia and reperfusion.
      ). However, the ATP/ADP ratio during reperfusion measured in guinea pig hearts could not be reproduced precisely (see Fig. 3F). Notably, in the work of Clarke et al. (
      • Clarke K.
      • O'Connor A.J.
      • Willis R.J.
      Temporal relation between energy metabolism and myocardial function during ischemia and reperfusion.
      ), this ratio was not directly measured but calculated using the creatine phosphate and pH measurements and assuming a constant creatine pool and creatine kinase equilibrium. Two possible reasons for the mismatch between model predictions and experimental results could be the following: (A) The difference in the total amount of creatine ([CrP] + [Cr]) results in different ATP/ADP dynamics during ischemia and reperfusion. Figure 3D shows that creatine phosphate is fully recovered after about 60 s according to the model, while the experimental values indicate that it only takes approx. 20 s. These different timescales can also be observed in Figure 3F. (B) The assumption of a constant cytosolic pH results in a too simplified model such that the ATP/ADP dynamics upon reperfusion cannot be captured correctly. However, in this context this is unlikely since Clarke et al. report that pH is the value with the slowest transition time during ischemia and reperfusion. It is also possible that the high collateral circulation present in guinea pig hearts results in a higher residual oxygen concentration than used in the simulations (Clarke et al. showed that there were no significant drops in ATP during 150 s of ischemia in guinea pig hearts). To verify this, simulations were also carried out with higher residual oxygen concentrations (up to 10 % of physiological levels); however, this did not result in an initial overshoot of the ATP/ADP ratio, as it can be seen in Figure 3F. Even though the characteristic overshoot of the ATP/ADP ratio seen in Figure 3F cannot be reproduced by the model, it accurately predicts that the ATP/ADP ratio reaches steady state after approx. 75 s. This is also consistent with the experimental data.
      The model predictions were also compared with experimental measurements by Wu et al. (
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      ), who used phosphate-magnetic resonance spectroscopy (P-MRS) to measure creatine phosphate concentrations in the left anterior descending artery in canine hearts under different cardiac work states. This comparison showed that predictions of creatine phosphate levels during ischemia and reperfusion correlate well with experimental values (see Fig. 2). Experimental and predicted data both show a close to instant response to reduction in cellular oxygen levels and use the creatine phosphate stores to keep ATP levels at preischemic values for as long as possible. Furthermore, these creatine phosphate stores are also refilled immediately, once sufficient oxygen is available to generate ATP through oxidative phosphorylation.
      Sensitivity analysis at multiple timepoints shows how the control mechanism of the metabolic network changes between physiological and pathological conditions. During normal oxygen concentrations, most of the glycolytic metabolite levels are regulated by the enzymes hexokinase and phosphofructokinase. This behavior is consistent with other studies: indeed, hexokinase, phosphofructokinase, and pyruvate kinase are generally assumed to be the flux controlling enzymes of glycolysis. However, it has been shown that pyruvate kinase does not seem to limit the overall flux of glycolysis in hearts since the rate of pyruvate kinase is linked to the glycolytic flux (
      • van Erp H.E.
      • Roholl P.J.
      • Rijksen G.
      • Sprengers E.D.
      • van Veelen C.W.
      • Staal G.E.
      Production and characterization of monoclonal antibodies against human type K pyruvate kinase.
      ,
      • Depré C.
      • Rider M.H.
      • Hue L.
      Mechanisms of control of heart glycolysis.
      ). Additionally, the concentrations of energy-related molecules such as adenosine and creatine are strongly regulated by the ATP consumption rate and the adenine nucleotide translocase (a transporter protein that exchanges ATP and ADP between cellular compartments).
      During ischemic conditions, mitochondrial reactions included in the model do not exhibit any control since oxidative phosphorylation needs oxygen to produce ATP. Therefore, all reaction rates in the mitochondria related to oxidative phosphorylation asymptotically converge to 0. Surprisingly, however, the enzyme exerting the highest regulation of cytosolic concentrations is glucose-6-phosphate dehydrogenase, which acts as a sink for the glycolysis intermediate G6P. This is interesting because the pentose phosphate pathway is not included in this model. The pentose phosphate pathway is a metabolic pathway parallel to glycolysis that produces sugars to make DNA, RNA, and NADPH molecules used in, e.g., fatty acid synthesis (
      • Wamelink M.M.
      • Struys E.A.
      • Jakobs C.
      The biochemistry, metabolism and inherited defects of the pentose phosphate pathway: A review.
      ).
      As expected, the sensitivity analysis performed after reperfusion shows similar results to preischemic conditions since the reaction parameters are not influenced by the ischemic period. This behavior arises due to the fact that there are no mechanisms for irreversibility or cellular damage included in the current model. An improved model version would potentially include a mechanism for cellular damage as a function of the time and severity of ischemia. It is common in ischemia/reperfusion simulations to change reaction parameters between pre- and postischemic periods (e.g., Chouchani et al. (
      • Chouchani E.T.
      • Pell V.R.
      • Gaude E.
      • Aksentijević D.
      • Sundier S.Y.
      • Robb E.L.
      • Logan A.
      • Nadtochiy S.M.
      • Ord E.N.J.
      • Smith A.C.
      • Eyassu F.
      • Shirley R.
      • Hu C.H.
      • Dare A.J.
      • James A.M.
      • et al.
      Ischaemic accumulation of succinate controls reperfusion injury through mitochondrial ROS.
      ) reduced the maximum reaction rate of ATP synthase to 50%, and Bazil et al. (
      • Bazil J.N.
      • Beard D.A.
      • Vinnakota K.C.
      Catalytic coupling of oxidative phosphorylation, ATP demand, and reactive oxygen species generation.
      ) increased the complex II activity and equilibrium constant by a factor of 10 after reperfusion). However, in many cases these parameter changes have only been qualitatively estimated, while in the current model, the focus was on determining the biochemical changes on a quantitative level.

      Oxygen dependence and oxidative damage

      Results shown in Figure 4 suggest that cardiomyocytes are able to survive at oxygen concentrations as low as 10% of the normal physiological value. This is in agreement with results reported by McDougal and Dewey (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ), who suggested a minimum viable oxygen concentration of 0.01 mM (8.5% of physiol. level) for nonbeating cardiomyocytes and 0.006 mM (4.6% of physiol. level) for cardiomyocytes that just keep up the most essential functions to survive. This is especially interesting, because McDougal and Dewey assumed a constant ATP consumption rate for the nonbeating and maintenance-level states of 56% and 21% of the healthy ATP consumption. In this work, we refined this assumption by including the dynamic ATP consumption rate suggested by Wu et al. (
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      ). Wu et al. modeled ATP consumption as a function of the relative molecular concentrations of Pi, ADP, and ATP and validated the predictions in full-occlusion protocols on canine hearts. Additionally, the values were consistent with similar protocols on canine hearts, measured by Schwartz et al. (
      • Schwartz G.G.
      • Schaefer S.
      • Meyerhoff D.J.
      • Gober J.
      • Fochler P.
      • Massie B.
      • Weiner M.W.
      Dynamic relation between myocardial contractility and energy metabolism during and following brief coronary occlusion in the pig.
      ). Therefore, the current model does not make any assumptions about the state of the cardiomyocyte (i.e., beating, nonbeating, maintenance-only), but implicitly adjusts the ATP consumption and consequently creates a potential readout for the cardiomyocyte state based on metabolic measurements. However, the parameters in the equation associated with the ATP consumption rate are based on experiments with canine hearts. Therefore, further biological experiments using human cardiomyocytes will be necessary to assess whether these results are also applicable to humans.
      During ischemia, the mitochondrial membrane potential asymptotically approaches a value of approximately 55 mV. Dissecting the equation for the change of the membrane potential over time into its components (see Table 2) shows that during prolonged ischemia, all chemical reaction rates approach 0 except for the proton leak. Thus, the membrane potential is only dissipated via the proton leak. This is also in agreement with Borutaite et al. (
      • Borutaite V.
      • Mildaziene V.
      • Brown G.C.
      • Brand M.D.
      Control and kinetic analysis of ischemia-damaged heart mitochondria: Which parts of the oxidative phosphorylation system are affected by ischemia?.
      ), who measured the effects of the respiratory chain and the proton leak on the mitochondrial membrane potential in isolated mitochondria. The authors showed that during ischemia, the control coefficient for the proton leak increases while the coefficient for the respiratory chain decreases. It was further hypothesized that this was driven by an increase in permeability of the mitochondrial membrane to protons.
      Table 2Rate of change for species in mitochondrial intermembrane space
      d[Species]dtVIMSpaceTypeExpression
      [ADP]iAssignment[ADP]x
      [fADP]iODEJMgADPi
      [mADP]iODEJMgADPi
      [ATP]iAssignment[ATP]x
      [fATP]iODEJMgATPi
      [mATP]iODEJMgATPi
      [H]iConst.n.A.
      [K]iConst.n.A.
      [Pi]iODEJPi2JPi1
      [Cox]iODE2JC42JC3
      [Cred]iODE2JC32JC4
      ΔΨODE(4JC1+2JC3+4JC43JF1JANTJHle)/CIM
      Except for ΔΨ, where the rate of change directly refers to dΔΨdt.
      Type indicates the method used to calculate concentration levels. ODE: expression refers to rate of change d[Species]dtVIMSpace. Assignment: concentration is directly calculated based on expression. Const: Concentration is fixed.
      Except for ΔΨ, where the rate of change directly refers to dΔΨdt.
      Although we have modeled only 15 min of ischemia, cardiomyocyte damage can occur during this short period. Early observations of cardiomyocytes revealed structural changes to the mitochondria after only 5 to 10 min of ischemia (
      • Schaper J.
      • Schaper W.
      Time course of myocardial necrosis.
      ), and more recent studies have shown depletion of desmin, actin, myoglobin, along with the presence of fibrinogen and C5 in cardiomyocytes after only 15 min of ischemia. Early canine studies also indicated that short periods of ischemia (15–20 min) were accompanied by substantial accumulation of lactic acid, along with depletion of glycogen and ATP (
      • Jennings R.B.
      Early phase of myocardial ischemic injury and infarction.
      ,
      • Schaper J.
      • Mulch J.
      • Winkler B.
      • Schaper W.
      Ultrastructural, functional, and biochemical criteria for estimation of reversibility of ischemic injury: A study on the effects of global ischemia on the isolated dog heart.
      ). While these effects could be reversed quickly during reperfusion, cardiac function was significantly impacted and slow to recover after these short ischemic intervals, as indicated by a high incidence of ventricular fibrillation upon reperfusion (
      • Jennings R.B.
      Early phase of myocardial ischemic injury and infarction.
      ), and persistently depressed contractile function even after 50 min of reperfusion (
      • Schaper J.
      • Mulch J.
      • Winkler B.
      • Schaper W.
      Ultrastructural, functional, and biochemical criteria for estimation of reversibility of ischemic injury: A study on the effects of global ischemia on the isolated dog heart.
      ). Furthermore, myocardial electrical impedance has been shown to increase in a biphasic manner after ischemia, exhibiting an early increase after onset of ischemia, a plateau phase, and a second rapid increase associated with cellular uncoupling (
      • del Rio C.L.
      • McConnell P.I.
      • Clymer B.D.
      • Dzwonczyk R.
      • Michler R.E.
      • Billman G.E.
      • Howie M.B.
      Early time course of myocardial electrical impedance during acute coronary artery occlusion in pigs, dogs, and humans.
      ). Myocardial electrical impedance reaches a plateau after 46 min of ischemia in dogs, but after less than 5 min of ischemia in humans, reflecting the observation that human myocardium is particularly sensitive to ischemia, even if it is of short duration (
      • del Rio C.L.
      • McConnell P.I.
      • Clymer B.D.
      • Dzwonczyk R.
      • Michler R.E.
      • Billman G.E.
      • Howie M.B.
      Early time course of myocardial electrical impedance during acute coronary artery occlusion in pigs, dogs, and humans.
      ).
      With decreasing ischemic O2 values, the ATP pool starts to deplete at an increasing speed (see Fig. 4A). Additionally, the data in Figure 4 indicates that, at O2 levels below 10% of normal physiological levels, conditions favorable for the generation of reactive oxygen species are present. Below this threshold value, the model predicts a reversal of the electron transport at complex I and a spike in the mitochondrial membrane potential, which is associated with excessive superoxide production at complex I and III, respectively (
      • Chouchani E.T.
      • Pell V.R.
      • Gaude E.
      • Aksentijević D.
      • Sundier S.Y.
      • Robb E.L.
      • Logan A.
      • Nadtochiy S.M.
      • Ord E.N.J.
      • Smith A.C.
      • Eyassu F.
      • Shirley R.
      • Hu C.H.
      • Dare A.J.
      • James A.M.
      • et al.
      Ischaemic accumulation of succinate controls reperfusion injury through mitochondrial ROS.
      ,
      • Padmaraj D.
      • Pande R.
      • Miller J.H.
      • Wosik J.
      • Zagozdzon-Wosik W.
      Mitochondrial membrane studies using impedance spectroscopy with parallel pH monitoring.
      ,
      • Korshunov S.S.
      • Skulachev V.P.
      • Starkov A.A.
      High protonic potential actuates a mechanism of production of reactive oxygen species in mitochondria.
      ). Both of these phenomena could potentially be explained by the altered redox state of the electron transport chain. In physiological conditions, the individual protein complexes in the electron transport chain act as a supercomplex and their activity is highly coupled. This is also observed in model predictions, where it can be seen that the reaction rates at complex I and III are highly correlated during the preischemic and ischemic period. However, this correlation is not observed anymore during reperfusion, as shown in Figure 4, B and D. During ischemia, the sequential stopping of complex IV, III, and then I leads to fully reduced coenzyme Q and cytochrome C pools. Upon reperfusion, the cytochrome C pool rapidly drives the proton pumping of complex IV leading to a spike in the membrane potential. Simultaneously, the redox driving force resulting from the reduced state of coenzyme Q leads to a reversal of the chemical reaction at complex I and therefore reverse electron transport (RET). The chemical flux through complex I is represented as a thermodynamically balanced model and is driven toward a steady state, which is a function of the proton gradient and the mitochondrial membrane potential. A reduced coenzyme Q pool during ischemia combined with a high mitochondrial membrane potential thus drives the chemical reaction at complex I to reverse and as a consequence to a reverse transfer of electrons from ubiquinol to complex I. This has also been discussed by Scialò et al. (
      • Scialò F.
      • Fernández-Ayala D.
      • Sanz A.
      Role of mitochondrial reverse electron transport in ROS signaling: Potential roles in health and disease.
      ), who reviewed different mechanisms of ROS production at complex I during forward and reverse electron transport.
      While the elevation of complex III activity (together with complex IV) is quite short-lived, the hyperpolarization during reperfusion negatively correlates with activity of the mitochondrial proton leak, suggesting that increased activity of the proton leak may impact the mitochondrial membrane potential. This has also been discussed by Cheng et al. (
      • Cheng J.
      • Nanayakkara G.
      • Shao Y.
      • Cueto R.
      • Wang L.
      • Yang W.Y.
      • Tian Y.
      • Wang H.
      • Yang X.
      Mitochondrial proton leak plays a critical role in pathogenesis of cardiovascular diseases.
      ), who mentioned the existence of a protective feedback loop between ROS generation and proton leak activity. Based on our model results, we hypothesize that this feedback loop may be involved in the membrane potential observed during reperfusion.
      Chouchani et al. recently proposed that RET is a key driver of oxidative damage during reperfusion (
      • Chouchani E.T.
      • Pell V.R.
      • Gaude E.
      • Aksentijević D.
      • Sundier S.Y.
      • Robb E.L.
      • Logan A.
      • Nadtochiy S.M.
      • Ord E.N.J.
      • Smith A.C.
      • Eyassu F.
      • Shirley R.
      • Hu C.H.
      • Dare A.J.
      • James A.M.
      • et al.
      Ischaemic accumulation of succinate controls reperfusion injury through mitochondrial ROS.
      ,
      • Chouchani E.T.
      • Pell V.R.
      • James A.M.
      • Work L.M.
      • Saeb-Parsy K.
      • Frezza C.
      • Krieg T.
      • Murphy M.P.
      A unifying mechanism for mitochondrial superoxide production during ischemia-reperfusion injury.
      ). Our model predicts a drastically increased mitochondrial membrane potential during the first minutes of reperfusion. A potential explanation for this could be the large redox driving force that pushes protons across the mitochondrial membrane to establish the electrochemical gradient used for ATP production. During ischemia, the coenzyme Q and cytochrome C pools become fully reduced. Thus, at the first moments of reperfusion, reactions at complex I, III, and IV will have a shifted chemical equilibrium resulting in faster reaction rates and accordingly more protons being pumped across the membrane simultaneously. It has been shown that an increased mitochondrial membrane potential can lead to a strong increase in superoxide production at complex III (
      • Padmaraj D.
      • Pande R.
      • Miller J.H.
      • Wosik J.
      • Zagozdzon-Wosik W.
      Mitochondrial membrane studies using impedance spectroscopy with parallel pH monitoring.
      ,
      • Korshunov S.S.
      • Skulachev V.P.
      • Starkov A.A.
      High protonic potential actuates a mechanism of production of reactive oxygen species in mitochondria.
      ). A high membrane potential slows down the electron transfer from heme bL to heme bH (
      • Bleier L.
      • Dröse S.
      Superoxide generation by complex III: From mechanistic rationales to functional consequences.
      ), two subunits of complex III. Therefore, the predicted spike of the membrane potential seems to lead to favorable conditions for superoxide production at complex III (
      • Rottenberg H.
      • Covian R.
      • Trumpower B.L.
      Membrane potential greatly enhances superoxide generation by the cytochrome bc1 complex reconstituted into phospholipid vesicles.
      ). In isolated mitochondria under state 4 conditions and succinate (supporting RET), mitochondrial membrane potential increases of 30 to 50 mV (to nearly 200 mV) have been observed (
      • Aon M.A.
      • Cortassa S.
      • O’Rourke B.
      Redox-optimized ROS balance: A unifying hypothesis.
      ). Although the absolute value of dPsi predicted is higher than that seen in biological experiments, the scenario of high dPsi, maximal NADH pool, and reverse electron transport has been previously described (
      • Aon M.A.
      • Cortassa S.
      • O’Rourke B.
      Redox-optimized ROS balance: A unifying hypothesis.
      ,
      • Miwa S.
      • Brand M.
      Mitochondrial matrix reactive oxygen species production is very sensitive to mild uncoupling.
      ). However, even though superoxide production due to the increase in membrane potential has been discussed in different studies (e.g., (
      • Suski J.
      • Lebiedzinska M.
      • Bonora M.
      • Pinton P.
      • Duszynski J.
      • Wieckowski M.R.
      Relation between mitochondrial membrane potential and ROS formation.
      )), few papers have been discussing the mechanism of superoxide production as a result of the mitochondrial membrane potential in the context of ischemia/reperfusion injury. It is therefore critical to design appropriate biological experiments that are able to confirm whether this identified mechanism is a driver of oxidative damage in the pathophysiology of reperfusion injury.

      A strategy to minimize ischemia/reperfusion injury

      Simulations of the stepwise reperfusion shown in Figures 5 and 6 suggest that biochemical modeling could help in identifying optimal reperfusion strategies to minimize tissue damage from reperfusion injury in clinical applications. Figure 5 indicates a correlation between tissue damage and severity of ischemia (i.e., recovery time for the mitochondrial membrane potential and time of reverse electron transport at complex I increase with lower oxygen concentrations during ischemia). Interestingly, however, the spiking of the membrane potential and the complex III reaction rate shown in Figure 4, C and D disappear when introducing the intermediate reperfusion step at 5% oxygen (see Fig. 5, B and D), suggesting an improved reperfusion outcome. In addition to the investigation of the relationship between tissue damage and severity of ischemia during the two-step reperfusion shown in Figures 5 and 6 indicates different degrees of oxidative damage as a function of the chosen reperfusion strategy. Thus, while Figure 5 is important from a biological perspective, Figure 6 is especially interesting from a clinical perspective as it directly suggests potential intervention strategies in a clinical setting. The mitochondrial membrane potential and complex I and III reaction rate shown in Figure 6 consistently indicate reduced damage from reperfusion injury when introducing an intermediate reperfusion step at optimized oxygen levels. The spike in the mitochondrial membrane potential (and the associated ROS production) is minimized when the oxygen concentration of the first reperfusion step is between 5% and 10% (see Fig. 6B, traces 3 and 4 from left to right). Similarly, Figure 6, C and D show a minimal RET at complex I and the spike in the reaction rate of complex III nearly disappeared when introducing this intermediate reperfusion step. These results are also supported by experimental findings of different groups (
      • Abdel-Rahman U.
      • Aybek T.
      • Moritz A.
      • Kleine P.
      • Matheis G.
      Graded reoxygenation limits lipid peroxidation during surgical reperfusion.
      ,
      • Douzinas E.E.
      • Pitaridis M.T.
      • Patsouris E.
      • Kollias S.
      • Boursinos V.
      • Karmpaliotis D.I.
      • Gratsias Y.
      • Evangelou E.
      • Papalois A.
      • Konstantinidou A.E.
      • Roussos C.
      Myocardial ischemia in intestinal postischemic shock: The effect of hypoxemic reperfusion.
      ,
      • Luo X.
      • Yin Y.
      • You G.
      • Chen G.
      • Wang Y.
      • Zhao J.
      • Wang B.
      • Zhao L.
      • Zhou H.
      Gradually increased oxygen administration improved oxygenation and mitigated oxidative stress after resuscitation from severe hemorrhagic shock.
      ). Abdel-Rahman et al. (
      • Abdel-Rahman U.
      • Risteski P.
      • Tizi K.
      • Kerscher S.
      • Behjati S.
      • Bejati S.
      • Zwicker K.
      • Scholz M.
      • Brandt U.
      • Moritz A.
      Hypoxic reoxygenation during initial reperfusion attenuates cardiac dysfunction and limits ischemia-reperfusion injury after cardioplegic arrest in a porcine model.
      ) showed in a porcine model and in a clinical study of 19 patients undergoing cardiac surgery (
      • Abdel-Rahman U.
      • Aybek T.
      • Moritz A.
      • Kleine P.
      • Matheis G.
      Graded reoxygenation limits lipid peroxidation during surgical reperfusion.
      ) that a graded reperfusion led to a decrease in myocardial oxidative injury. In a recent review of hypoxemic reperfusion, Tasoulis et al. (
      • Tasoulis M.K.
      • Douzinas E.E.
      Hypoxemic reperfusion of ischemic states: An alternative approach for the attenuation of oxidative stress mediated reperfusion injury.
      ) mention that a stepwise reperfusion reduces the available oxygen for ROS production, while still providing enough oxygen for the cell to recover from ischemia. Thus, the initial burst in ROS generation could be potentially mitigated by a gradual reperfusion. This indicates that (a) a multistep reperfusion strategy could outperform a one-step reperfusion regardless of ischemic oxygen concentrations and (b) a mathematical model of the cardiomyocyte metabolism, such as the one developed herein, could be used to identify optimal reperfusion strategies for clinical settings. To the best of our knowledge, the changes in molecular-level species during stepwise reperfusion and the associated biological mechanisms have not been studied extensively in vivo yet. Thus, it will be critical to design and perform adequate experiments that validate these findings and additionally to come up with an according theory that is able to explain the mechanisms responsible for the potentially beneficial effect of a stepwise reperfusion protocol.
      For the future, it will be important to refine the model in these directions to quantify the exact amounts of ROS being produced at the protein complexes of the electron transport chain and follow up these predictions with biological experiments. Gauthier et al. (
      • Gauthier L.D.
      • Greenstein J.L.
      • Cortassa S.
      • O’Rourke B.
      • Winslow R.L.
      A computational model of reactive oxygen species and redox balance in cardiac mitochondria.
      ) created such a computational model of ROS dynamics in cardiac mitochondria and suggest that ROS production increases exponentially as a function of the mitochondrial membrane potential. This further highlights the importance of our hypothesis that the membrane potential seems to be a key mechanism in ischemia/reperfusion injury.
      It is important to note that, under certain circumstances, ROS can also exert cardioprotective mechanisms by activating HIF-1 (
      • Cadenas S.
      ROS and redox signaling in myocardial ischemia-reperfusion injury and cardioprotection.
      ,
      • Eckle T.
      • Köhler D.
      • Lehmann R.
      • El Kasmi K.
      • Eltzschig H.K.
      Hypoxia-inducible factor-1 is central to cardioprotection: A new paradigm for ischemic preconditioning.
      ) and Nrf2 (
      • Zhang Y.
      • Sano M.
      • Shinmura K.
      • Tamaki K.
      • Katsumata Y.
      • Matsuhashi T.
      • Morizane S.
      • Ito H.
      • Hishiki T.
      • Endo J.
      • Zhou H.
      • Yuasa S.
      • Kaneda R.
      • Suematsu M.
      • Fukuda K.
      4-hydroxy-2-nonenal protects against cardiac ischemia-reperfusion injury via the Nrf2-dependent pathway.
      ) signaling pathways. In fact, the cardioprotective effects of ischemic preconditioning are mediated in part by activation of HIF-1, which inhibits mPTP opening (
      • Cai Z.
      • Zhong H.
      • Bosch-Marce M.
      • Fox-Talbot K.
      • Wang L.
      • Wei C.
      • Trush M.A.
      • Semenza G.L.
      Complete loss of ischaemic preconditioning-induced cardioprotection in mice with partial deficiency of HIF-1 alpha.
      ), and Nrf2 signaling, which increases activation and expression of antioxidant enzymes (
      • Angeloni C.
      • Motori E.
      • Fabbri D.
      • Malaguti M.
      • Leoncini E.
      • Lorenzini A.
      • Hrelia S.
      H2O2 preconditioning modulates phase II enzymes through p38 MAPK and PI3K/Akt activation.
      ,
      • Xu B.
      • Zhang J.
      • Strom J.
      • Lee S.
      • Chen Q.M.
      Myocardial ischemic reperfusion induces de novo Nrf2 protein translation.
      ). Additionally, postconditioning approaches, which bear some similarity to the two-step reperfusion protocol that we implement in our model, elicit cardioprotection through various mechanism that include attenuation of oxidative stress and inhibition of mPTP opening (
      • Hausenloy D.
      • Yellon D.
      Ischaemic conditioning and reperfusion injury.
      ,
      • Cohen M.
      • Yang X.
      • Downey J.
      The pH hypothesis of postconditioning: Staccato reperfusion reintroduces oxygen and perpetuates myocardial acidosis.
      ,
      • Argaud L.
      • Gateau-Roesch O.
      • Raisky O.
      • Loufouat J.
      • Robert D.
      • Ovize M.
      Postconditioning inhibits mitochondrial permeability transition.
      ,
      • Hausenloy D.
      • Ong S.
      • Yellon D.
      The mitochondrial permeability transition pore as a target for preconditioning and postconditioning.
      ,
      • Heusch G.
      • Boengler K.
      • Schulz R.
      Inhibition of mitochondrial permeability transition pore opening: The holy grail of cardioprotection.
      ). Although many studies have indicated that postconditioning reduces ROS production at reperfusion (
      • Kin H.
      • Zhao Z.
      • Sun H.
      • Wang N.
      • Corvera J.
      • Halkos M.
      • Kerendi F.
      • Guyton R.
      • Vinten-Johansen J.
      Postconditioning attenuates myocardial ischemia-reperfusion injury by inhibiting events in the early minutes of reperfusion.
      ,
      • Sun H.-Y.
      • Wang N.-P.
      • Kerendi F.
      • Halkos M.
      • Kin H.
      • Guyton R.A.
      • Vinten-Johansen J.
      • Zhao Z.-Q.
      Hypoxic postconditioning reduces cardiomyocyte loss by inhibiting ros generation and intracellular ca2+ overload.
      ), others have demonstrated that ROS scavenging attenuates the effects of postconditioning (
      • Tsutsumi Y.
      • Yokoyama T.
      • Horikawa Y.
      • Patel H.
      Reactive oxygen species trigger ischemic and pharmacological postconditioning: In vivo and in vitro characterization.
      ,
      • Penna C.
      • Rastaldo R.
      • Mancardi D.
      • Raimondo S.
      • Cappello S.
      • Gattullo D.
      • Losano G.
      • Pagliaro P.
      Postconditioning induced cardioprotection requires signaling through a redox-sensitive mechanism, mitochondrial ATP-sensitive K+ channel and protein kinase c activation.
      ,
      • Babiker F.
      • Al-Jarallah A.
      • Joseph S.
      Understanding pacing postconditioning-mediated cardiac protection: A role of oxidative stress and a synergistic effect of adenosine.
      ). Whether ROS acts as a mediator of injury or protection seems to depend critically on the amount and biological context of ROS production, with high or increased levels being associated with deleterious (
      • Cadenas S.
      ROS and redox signaling in myocardial ischemia-reperfusion injury and cardioprotection.
      ,
      • Chouchani E.T.
      • Pell V.R.
      • James A.M.
      • Work L.M.
      • Saeb-Parsy K.
      • Frezza C.
      • Krieg T.
      • Murphy M.P.
      A unifying mechanism for mitochondrial superoxide production during ischemia-reperfusion injury.
      ,
      • Muntean D.M.
      • Sturza A.
      • Dănilă M.D.
      • Borza C.
      • Duicu O.
      • Mornoş C.
      The role of mitochondrial reactive oxygen species in Cardiovascular injury and protective strategies.
      ) and low or minimal levels being associated with protective (
      • Tullio F.
      • Angotti C.
      • Perrelli M.-G.
      • Penna C.
      • Pagliaro P.
      Redox balance and cardioprotection.
      ) effects. Although we do not directly model ROS generation in our system, we believe that the large change in dPsi during reperfusion would generate relatively high levels of ROS, since ROS production at complex III has exponential dependence on dPsi (
      • Rottenberg H.
      • Covian R.
      • Trumpower B.L.
      Membrane potential greatly enhances superoxide generation by the cytochrome bc1 complex reconstituted into phospholipid vesicles.
      ,
      • Guillaud F.
      • Dröse S.
      • Kowald A.
      • Brandt U.
      • Klipp E.
      Superoxide production by cytochrome bc1 complex: A mathematical model.
      ). We therefore hypothesize that the resultant high levels of ROS would exert a deleterious, rather than protective, effect on the myocyte.
      The model could be extended to include other energy-related pathways such as fatty acid oxidation and processes related to muscular contractions including calcium dynamics and more general electrophysiological mechanisms (
      • Winslow R.L.
      • Cortassa S.
      • O'Rourke B.
      • Hashambhoy Y.L.
      • Rice J.J.
      • Greenstein J.L.
      Integrative modeling of the cardiac ventricular myocyte.
      ,
      • Cortassa S.
      • Aon M.A.
      • Sollott S.J.
      Control and regulation of substrate selection in cytoplasmic and mitochondrial catabolic networks. A systems biology analysis.
      ). With our model being primarily focused on glucose as a substrate, this could be an excellent addition to study the role of additional metabolic substrates within the scope of ischemia/reperfusion injury.
      Yaniv et al. (
      • Yaniv Y.
      • Juhaszova M.
      • Nuss H.B.
      • Wang S.
      • Zorov D.B.
      • Lakatta E.G.
      • Sollott S.J.
      Matching ATP supply and demand in mammalian heart: In vivo, in vitro, and in silico perspectives.
      ) highlighted ADP/Pi and calcium as main regulators of ATP production and consumption in the mammalian heart. While the ADP and Pi controlled ATP consumption rate was included in our work, we do not have a direct relationship between calcium and ATP consumption. Calcium at least partially controls mitochondrial ATP production by regulating the activation of various enzymes in the Krebs cycle and the electrochemical gradient established in the electron transport chain. In order to capture these mechanisms, the phenomenological representation of the Krebs cycle used in this model could be replaced by including the chemical kinetics of the individual dehydrogenases, such as done by Wu et al. (
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      ) in their mathematical representation of the Krebs cycle.
      Our model does not yet incorporate mechanisms that directly tie changes in mitochondrial injury during reperfusion to cell death mechanisms. It has been shown that excessive ROS production, rapid pH restoration, and intracellular calcium overload during reperfusion (
      • Hausenloy D.J.
      • Yellon D.M.
      Myocardial ischemia-reperfusion injury: A neglected therapeutic target.
      ) can lead to disruption of mitochondrial integrity through the opening of mitochondrial permeability transition pore (MPTP) (
      • Gateau-Roesch O.
      • Argaud L.
      • Ovize M.
      Mitochondrial permeability transition pore and postconditioning.
      ,
      • Lesnefsky E.J.
      • Chen Q.
      • Tandler B.
      • Hoppel C.L.
      Mitochondrial dysfunction and myocardial ischemia-reperfusion: Implications for novel therapies.
      ). MPTP opening, in turn, can cause mitochondrial swelling and rupture, releasing intermembrane proteins such as cytochrome c and apoptosis-inducing factor that activate cell death pathways (
      • Xu A.
      • Szczepanek K.
      • Hu Y.
      • Lesnefsky E.J.
      • Chen Q.
      Cardioprotection by modulation of mitochondrial respiration during ischemia-reperfusion: Role of apoptosis-inducing factor.
      ). Several types of cell death, including apoptosis and necrosis, are implicated in I/R injury (
      • Hotchkiss R.S.
      • Strasser A.
      • McDunn J.E.
      • Swanson P.E.
      Cell death.
      ,
      • Del Re D.P.
      • Amgalan D.
      • Linkermann A.
      • Liu Q.
      • Kitsis R.N.
      Fundamental mechanisms of regulated cell death and implications for heart disease.
      ), and future work should focus on expanding our computational framework to model these mechanisms.
      In the scope of this work, only the temporal dynamics of the cardiac metabolism were incorporated. However, the concentration of many species and organelles such as the mitochondria are not spatially homogeneous. Balaban et al. (
      • Glancy B.
      • Hartnell L.M.
      • Combs C.A.
      • Femnou A.
      • Sun J.
      • Murphy E.
      • Subramaniam S.
      • Balaban R.S.
      Power grid protection of the muscle mitochondrial reticulum.
      ) showed that cardiac mitochondria are organized in regional subnetworks, which was suggested to limit the spread of mitochondrial dysfunction to local regions of the cell. In addition, the model could be extended to include the time delays inherent in reperfusion dynamics and diffusion of oxygen from the vasculature to the cardiomyocytes, as both effects would influence the time course of oxygen delivery to the reperfused cardiac tissues. To include such mechanisms, an improved simulation of the cardiac metabolism would be based on a spatiotemporal model.
      However, the focus of this work was to develop an experimentally validated model of the cardiomyocyte metabolism to gain insights into the complex processes and reactions present during ischemia/reperfusion injury. Additionally, we aim to provide the research community with a metabolic model that can be easily extended based on a specific research question, such as ROS production, cardiac muscle contractions, or fatty acid oxidation. To achieve this, the herein developed model is available online on http://www.ebi.ac.uk/biomodels/. Modeling cellular conditions though ischemia and reperfusion will enable researchers to study the mechanisms of reperfusion injury in cardiomyocytes in silico and also help to plan and conduct biological experiments in order to improve the understanding of the pathophysiology of ischemia/reperfusion injury.
      In summary, we have implemented a model of the cardiomyocyte metabolism, which is able to simulate the pathophysiological states of ischemia/reperfusion injury as a function of different extracellular oxygen levels. Additionally, the model was qualitatively and quantitatively validated using experimental measurements several animal species. The model suggests that cardiomyocytes are capable of adjusting their energy production in response to drastically varying extracellular oxygen availabilities. More specifically, the results show that oxygen concentrations as low as 10% of the normal physiological levels are enough for the cells to survive. Our results suggest that changes in the redox state of the electron transport chain are likely responsible for oxidative damage occurring during reperfusion injury. Reverse electron transport and a highly increased mitochondrial membrane potential both create conditions for increased rates of ROS production and thus oxidative damage. Most importantly, the model shows potential intervention strategies such as a two-step reperfusion protocol, in order to minimize oxidative damage during the first moments of reperfusion.

      Experimental procedures

      To track the time course of ATP and other related species, a model comprised of five compartments and 67 species was created using a systems biology approach. This resulted in a nonlinear system of ordinary differential equations (ODEs). This implementation builds upon the model McDougal and Dewey (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ), describing the anaerobic metabolism in cardiomyocytes. Their model is mostly focused on the cytosolic part of the energy production, i.e., glycolysis, glycogenolysis, and glycogenesis and cytosolic ATP buffering. In this work, this description was extended by incorporating a more complete and exact model of the mitochondrial ATP production inspired by the work of Beard (
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ). Additionally, the constant ATP consumption rate was replaced by a phosphate and ADP-dependent rate as suggested by Wu et al. (
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      ). Generally speaking, the model starts at the diffusion of oxygen from the blood vessel into the extracellular space. From there, glucose and oxygen diffuse through the cellular membrane into the cytosol of the cardiomyocyte. Glucose then enters the glycolytic pathway with pyruvate being the final metabolite. In our model, we integrated the conversion of pyruvate into acetyl-coenzyme A and the subsequent citric acid cycle into a single phenomenological process that reduces NAD to NADH, as suggested by Beard (
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ). This NADH then enters the electron transport chain at complex I to transfer electrons to complex IV with oxygen being the final electron acceptor. Simultaneously these protein complexes transfer protons from the mitochondrial matrix into the intermembrane space. The resulting proton-motive force, consisting of the proton gradient and the mitochondrial membrane potential, is finally used by the F1F0-ATP synthase to produce ATP from ADP and inorganic phosphate.
      To formulate a system of ODEs representing the dynamics of the metabolic processes, the reactions involving this species were listed for every species represented in the model. Afterward, the flux of each reaction and/or diffusion process was calculated based on either Michaelis–Menten kinetics or diffusion equations. Finally, the rate of change for each species was calculated as the sum of the reactions producing this species minus the sum of reactions consuming the species plus (or minus) the reactions/diffusion processes that transport a species in or out of its compartment and represented by a system of nonlinear ordinary differential equations. This is illustrated in the following equation
      d[Ci]dt=JProduction±JTransportJConsumption


      where [Ci] refers to the concentration of species i in mM, JProduction and JConsumption represent the amount of [Ci] that is being produced and consumed respectively and JTransport is the amount of species i that is transported in and out of the compartment per unit time. For a full list of the model equations and initial values, see the corresponding tables in the following section.

      Compartments

      Cytosol

      Most species in the cytosol are: (1) associated to glycolysis; (2) forms of adenosine phosphate; (3) related buffer systems. All kinetic parameters related to processes in the cytosol are taken from McDougal and Dewey (
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ), except transport processes between the cytosol and the mitochondria. See Table 3 for species and initial conditions and Table 4 for the rate of change for each species.
      Table 3Cytosolic species
      SpeciesFull nameUnitInitial valueSource
      [AMP]cAdenosine MonophosphatemM1×105(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [ADP]cAdenosine DiphosphatemM1×105(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [ATP]cAdenosine TriphosphatemM7(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [NAD]cNicotinamide Adenine Dinucleotide (ox)mM2.2565(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [NADH]cNicotinamide Adenine Dinucleotide (red)mM0.7135(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [Glucose]cGlucosemM1.91(
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      )
      [13BPG]c1,3-BisphosphoglyceratemM8.69×104(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [2PG]c2-PhosphoglyceratemM0.009(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [3PG]c3-PhosphoglyceratemM0.071(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [Cr]cCreatinemM0(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [CrP]cCreatine PhosphatemM25(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [DHAP]cDihydroxyacetone PhosphatemM0.036(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [F16BP]cFructose 1,6-BisphosphatemM6.78×104(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [F6P]cFructose-6-PhosphatemM0.041(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [G16BP]cGlucose-1,6-BisphosphatemM0.007(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [G1P]cGlucose-1-PhosphatemM0.02(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [G6P]cGlucose-6-PhosphatemM0.169(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [GAP]cGlyceraldehyde-3-PhosphatemM0.00162(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [Glycogen]cGlycogenmM21.4(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [Li]cLactatemM0.247(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [LA]cLactic AcidmM155.84(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [Mb]cMyoglobinmM0.00543(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [MbO2]cSaturated MyoglobinmM0.18457(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [O2]cOxygenmM0.11017(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [PEP]cPhosphoenolpyruvatemM0.013(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [pH]cpHmM7.1(
      US Burden of Disease Collaborators
      The state of US health, 1990-2010: Burden of diseases, injuries, and risk factors..
      )
      [Pi]cInorganic PhosphatemM7(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      [PYR]cPyruvatemM0.055(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [UDPG]cUridine Diphosphate GlucosemM0.099(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      [H]cHydrogenmM7.94×104(
      US Burden of Disease Collaborators
      The state of US health, 1990-2010: Burden of diseases, injuries, and risk factors..
      )
      [K]cPotassiummM150(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Mg]cMagnesiummM5(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [mADP]cMagnesium-bound ADPmM9.35×106(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      Table 4Rate of change for cytosolic species
      d[Species]dtVcytosolTypeExpression
      [AMP]cODEJAK
      [ADP]cODEJCK2JAK+JHK+JPFKJPGKJPK+JATPaseαF1JF1
      [ATP]cODEJCK+JAKJATPaseJHKJPFK+JPGK+JPK+αF1JF1
      [NAD]cConst.n.A.
      [NADH]cConst.n.A.
      [Glucose]cConst.n.A.
      [13BPG]cODEJGAPDHJPGKJG16BPS
      [2PG]cODEJPGMJEnolase
      [3PG]cODEJPGKJPGM+JG16BPS
      [Cr]cODEJCK
      [CrP]cODEJCK
      [DHAP]cODEJFBPAJTPIJG3PDH
      [F16BP]cODEJPFKJFBPA
      [F6P]cODEJPGIJPFK
      [G16BP]cConst.n.A.
      [G1P]cODEJPGluMJUDPGP+JGPJG16BPS
      [G6P]cODEJHKJPGI+JPGluMJG6PDH
      [GAP]cODEJFBPA+JTPIJGAPDH
      [Glycogen]cODEJGSD+JGSIJGP
      [Li]cODEJLFiJLFoJLA+JLDH
      [LA]cODEJLA
      [Mb]cODEJO2Mb
      [MbO2]cODEJO2Mb
      [O2]cODEJC4+JO2,ECJO2Mb
      [PEP]cODEJEnolaseJPK
      [Pi]cODEJGAPDH+2JUDPGPJGP+JG6PDH+JG3PDH+2JG16BPS+JATPaseαF1JF1
      [PYR]cODEJPKJLDHαDHJDH
      [UDPG]cODEJUDPGPJGSDJGSI
      [H]cConst.n.A.
      [K]cConst.n.A.
      [Mg]cAssignment[Mg]c,total[mADP]c
      [mADP]cAssignmentβADP
      Type indicates the method used to calculate concentration levels. ODE: expression refers to rate of change d[Species]dtVcytosol. Assignment: concentration is directly calculated based on expression. Const: Concentration is fixed.

      Mitochondrial intermembrane space

      The intermembrane space acts mostly as a storage to sustain the proton gradient and the mitochondrial membrane potential necessary for the F1F0-ATP synthase to produce ATP from ADP and inorganic phosphate. Furthermore, since the F1F0-ATP synthase requires ADP and magnesium as a cofactor, binding between magnesium ions and ATP and ADP is also included in the intermembrane space and the mitochondrial matrix. Here, ADP and ATP always refer to the sum of magnesium bound (mADP, mATP) and free species (fADP, fATP). See Tables 2 and 5 for all species and corresponding rate of change equations in the intermembrane space.
      Table 5Species in mitochondrial intermembrane space
      SpeciesFull nameUnitInitial valueSource
      [ADP]iAdenosine DiphosphatemM10(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [fADP]ifree ADPmM10(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [mADP]iMagnesium-bound ADPmM0(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [ATP]iAdenosine TriphosphatemM7(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [fATP]ifree ATPmM7(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [mATP]iMagnesium-bound ATPmM0(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [H]iHydrogenmM7.94×105(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [K]iPotassiummM150(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Pi]iInorganic PhosphatemM1(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Cox]iCytochrome C (ox)mM1.7(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Cred]iCytochrome C (red)mM1(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      ΔΨMitochondrial Membrane PotentialmV170(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )

      Mitochondrial matrix

      All processes in the mitochondrial matrix are either related to oxidative phosphorylation or transport processes between the mitochondrial matrix and the intermembrane space. First, pyruvate enters the mitochondrial matrix. Here it is worth mentioning that in this model the conversion of pyruvate into acetyl-coenzyme A and the TCA cycle are modeled as a single phenomenological dehydrogenase flux converting NAD into NADH, as suggested by Beard (
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ). This NADH is then used by the first protein complex (CI) in the electron transport chain to pump ions across the mitochondrial membrane. Similarly complex III and IV are pumping ions across the membrane to establish a proton-motive force, which is finally used by the F1F0-ATP synthase to produce ATP from ADP and inorganic phosphate. Although we have not modeled complex II explicitly, we incorporate the behavior of this complex implicitly via reactions occurring at complex I and III via changes in the Q and QH2 pools. For an overview of the species and rate of change equations in the mitochondrial matrix, refer to Tables 6 and 7, respectively.
      Table 6Species in mitochondrial matrix
      SpeciesFull nameUnitInitial valueSource
      [ADP]mAdenosine DiphosphatemM10(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [fADP]mfree ADPmM10(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [mADP]mMagnesium-bound ADPmM0(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [ATP]mAdenosine TriphosphatemM7(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [fATP]mfree ATPmM7(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [mATP]mMagnesium-bound ATPmM0(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [H]mHydrogenmM6.31×105(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [K]mPotassiummM140(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Mg]mMagnesiummM5(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [NAD]mNicotinamide Adenine Dinucleotide (ox)mM1.47(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [NADH]mNicotinamide Adenine Dinucleotide (red)mM1.5(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Pi]mInorganic PhosphatemM1(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [Q]mUbiquinonemM0.55(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      [QH2]mUbiquinolmM0.8(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      Table 7Rate of change for species in mitochondrial matrix
      d[Species]dtVMitochondriaTypeExpression
      [ADP]mODEJANTJF1
      [fADP]mODEJMgADPx
      [mADP]mODEJMgADPx
      [ATP]mODEJF1JANT
      [fATP]mODEJMgATPx
      [mATP]mODEJMgATPx
      [H]mODE1rbuff[H]m(JDH5JC12JC34JC4+2JF1+2JPi1+JHleJKH)
      [K]mODEJKH
      [Mg]mODEJMgADPxJMgATPx
      [NAD]mODEJC1JDH
      [NADH]mODEJDHJC1
      [Pi]mODEJPi1JF1
      [Q]mODEJC3JC1
      [QH2]mODEJC1JC3

      Extracellular space and blood vessel

      Species present in the extracellular space and the blood vessel and their rates of change are listed in Tables 8 and 9. Additionally, all of the model compartments and their relative sizes are provided in Table 10.
      Table 8Species in extracellular space and vessel
      SpeciesCompartmentFull nameUnitInitial valueSource
      [O2]eExtracellular SpaceOxygenmM0.1289(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [O2]vVesselOxygenmM0.1326(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      [L]eExtracellular SpaceLactatemM0.33(
      US Burden of Disease Collaborators
      The state of US health, 1990-2010: Burden of diseases, injuries, and risk factors..
      )
      Table 9Species in extracellular space and vessel
      d[Species]dtTypeExpression
      [O2]eODEJO2,VEJO2,ECVextraceullular
      [O2]vConst.n.A.
      Table 10Compartments involved in metabolic model
      CompartmentSizeUnitSource
      Vessel0.06842l(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Extracellular Space0.24063l(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Cytosol1l(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Mitochondrial Intermembrane Space0.0715l(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      Mitochondrial Matrix0.6435l(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )

      Simulations

      Time course

      Time course simulations were generally performed for the duration of 3000 s, where the first 1000 s corresponded to normal physiological conditions of the cell, the time between 1000 s and 2000 s represented ischemia by decreasing the vascular oxygen levels from 0.133 mM to values between 10 and 0.1% of the physiological levels and after 2000 s reperfusion was simulated by switching the oxygen levels in the blood vessel back to preischemic values.

      Sensitivity analysis

      To analyze the sensitivity of the concentration of each species with respect to the rate of each individual chemical reaction or transport process, the sensitivity coefficients were computed as follows:
      sij=d[Ci]dXj


      where [Ci] is the concentration of species i and Xj is the activity of the enzyme in reaction j. Furthermore, to work with unitless sensitivity coefficients, which allow for easier comparison between variables of different units (e.g., membrane potential and ATP concentration), the relative sensitivity coefficients s˜ij were calculated as
      s˜ij=sijXj[Ci]=d[Ci]dXjXj[Ci]


      In order to characterize the total influence of an enzyme on the molecular concentrations, the L1-norm of the relative sensitivities across all species was calculated as
      S˜j=iSpecies|s˜ij|


      Equations

      In the remainder of this section, an exhaustive list of chemical reactions and equations for reaction rates used in the model is given. The data are segmented into reactions present in the cytosol, intermembrane space, mitochondrial matrix, extracellular matrix, and the blood vessels surrounding the cell.

      Cytosolic enzymes and reactions

      ATP buffering

      ADP+CrP+HATP+Cr
      (R1)


      JCK=XCK(KCK[ADP]c[CrP]c[H]c[ATP]c[Cr]c)
      (1)


      Tabled 1
      ParameterValueUnitSource
      XCK10,000mM/(smM2)(
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      )
      KCK1,660,0001/mM(
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      )
      2ADPAMP+ATP
      (R2)


      JAK=XAK(KAK[ADP]c2[ATP]c[AMP]c)
      (2)


      Tabled 1
      ParameterValueUnitSource
      XAK10,000mM/(smM2)(
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      )
      KAK1unitless(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )

      Energy consumption

      ATPADP+Pi
      (R3)


      JATPase=XATPase1+R[Pi]c[ADP]c[ATP]c
      (3)


      Tabled 1
      ParameterValueUnitSource
      XATPase0.39mM/s(
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      )
      R65.81/mM(
      • Wu F.
      • Zhang E.Y.
      • Zhang J.
      • Bache R.J.
      • Beard D.A.
      Phosphate metabolite concentrations and ATP hydrolysis potential in normal and ischaemic hearts.
      )

      Glycolysis

      Glucose+ATPG6P+ADP+H
      (R4)


      JHK=Vf[Glucose]c[ATP]cKm,f([ATP]c+Km,ATP)Vr[G6P]cKm,r1+[Glucose]cKm,f+[G6P]cKm,r
      (4)


      Tabled 1
      ParameterValueUnitSource
      Vf0.550mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.072mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,ATP0.236mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr1.06×104mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.042mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      G6PF6P
      (R5)


      JPGI=Vf[G6P]cKm,fVr[F6P]cKm,r1+[G6P]cKm,f+[F6P]cKm,r
      (5)


      Tabled 1
      ParameterValueUnitSource
      Vf10.067mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.425mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr9.6mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.175mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      F6+ ATPF16+ ADP + H
      (R6)


      JPFK=Vf[F6P]cKm,f+[F6P]c1+Km,ATP[ATP]c
      (6)


      Tabled 1
      ParameterValueUnitSource
      Vf1.328mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.224mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,ATP0.127mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      F16+ ATPDHAP+GAP
      (R7)


      JFBPA=Vf[F16BP]c[F16BP]c+Km,f
      (7)


      Tabled 1
      ParameterValueUnitSource
      Vf0.992mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.038mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      DHAPGAP
      (R8)


      JTPI=Vf[DHAP]c[DHAP]c+Km,f
      (8)


      Tabled 1
      ParameterValueUnitSource
      Vf5.933mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f1.53mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      GAP+NAD+Pi13BPG+NADH+H
      (R9)


      JGAPDH=Vf[GAP]cKm,f+[GAP]c1+Km,NAD[NAD]c
      (9)


      Tabled 1
      ParameterValueUnitSource
      Vf5.35mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.042mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,NAD0.058mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      13BPG + ADP + H3PG + ATP
      (R10)


      JPGK=Vf[13BPG]cs1([13BPG]c)Km,f(1+Km,ADP[ADP]c)Vr[3PG]cs1([3PG]c)Km,r(1+Km,ATP[ATP]c)1+[13BPG]cKm,f+[3PG]cKm,r
      (10)


      Tabled 1
      ParameterValueUnitSource
      Vf251mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.021mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,ADP0.565mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr15.98mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.51mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,ATP0.008mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      3PG2PG
      (R11)


      JPGM=Vf[3PG]cs1([3PG]c)Km,fVr[2PG]cs1([2PG]c)Km,r1+[3PG]cKm,f+[2PG]cKm,r
      (11)


      Tabled 1
      ParameterValueUnitSource
      Vf11.233mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.145mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr48.0mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.139mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      2PGPEP+H2O
      (R12)


      JEnolase=Vf[2PG]s1([2PG])Km,fVr[PEP]Km,r1+[2PG]Km,f+[PEP]Km,r
      (12)


      Tabled 1
      ParameterValueUnitSource
      Vf1.85mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.045mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr2.00mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.089mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      PEP+ADPPYR+ATP
      (R13)


      JPK=Vf[PEP]cs1([PEP]c)Km,f(1+Km,ADP[ADP]c)Vr[PYR]cs1([PYR]c)Km,r1+[PEP]cKm,f+[PYR]cKm,r
      (13)


      Tabled 1
      ParameterValueUnitSource
      Vf9.433mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.11mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,ADP0.00268mM(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      ,
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr0.00105mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r10mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      G6Pø
      (R14)


      JG6PDH=Vf1vss,PGIJPGI
      (14)


      Tabled 1
      ParameterValueUnitSource
      Vf0.095mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      vss,PGI0.125mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      DHAP+NADHGlycerol3P+NAD
      (R15)


      JG3PDH=Vf[DHAP]c[NADH]cKia,NADHKm,DHAP+Km,DHAP[NADH]c+Km,NADH[DHAP]c+[DHAP]c[NADH]c
      (15)


      Tabled 1
      ParameterValueUnitSource
      Vf0.095mM/s(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Kia,NADH0.095mM(
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: A phantom phenomenon.
      )
      Km,DHAP0.095mM(
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: A phantom phenomenon.
      )
      Km,NADH0.095mM(
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: A phantom phenomenon.
      )
      13BPG+G1P3PG+G16BP
      (R16)


      JG16BPS=Vf[13BPG]cKm,f(1+Km,G1P[G1P]c)Vr[3PG]cKm,r(1+Km,G16BP[G16BP]c)1+[13BPG]cKm,f+[3PG]cKm,r
      (16)


      Tabled 1
      ParameterValueUnitSource
      Vf10mM/s(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Km,f0.021mM(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Km,G1P0.008mM(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Vr6mM/s(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Km,r0.51mM(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )
      Km,G16BP0.565mM(
      • McDougal A.D.
      • Dewey Jr., C.F.
      Modeling oxygen requirements in ischemic cardiomyocytes.
      )

      Glycogenesis and glycogenolysis

      G1PG6P
      (R17)


      JPGluM=Vf[G1P]cKm,fVr[G6P]cKm,r1+[G1P]cKm,f+[G6P]cKm,r
      (17)


      Tabled 1
      ParameterValueUnitSource
      Vf1.933mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.045mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr1.12mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r0.67mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      G1+ UTP + HUDPG+2Pi
      (R18)


      JUDPGP=XUDPGPkf[G1P]c[UDPG]ckr1+(Km,AMP[AMP]c)1.5
      (18)


      Tabled 1
      ParameterValueUnitSource
      XUDPGP10,000unitless(
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      )
      kf4.361/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,AMP0.016mM(
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      )
      kr0.8811/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      UDPGUDP+Glycogen
      (R19)


      JGSD=Vf,GSD[UDPG]c(Km,f,GSD+[UDPG]c)(1+(Km,AMP[AMP]c)1.5)
      (19)


      JGSI=Vf,GSI[UDPG]c(Km,f,GSI+[UDPG]c)(1+(Km,AMP[AMP]c)1.5)
      (20)


      Tabled 1
      ParameterValueUnitSource
      Vf,GSD0.147mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f,GSD1.42mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vf,GSI0.147mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f,GSI0.08mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,AMP0.016mM(
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      )
      Glycogen+PiG1P
      (R20)


      JGlyPhos=Vf[Glycogen]cKm,fVr[G1P]cKm,r(1+[Glycogen]cKm,f+[G1P]cKm,r)(1+(Km,AMP[AMP])1.5)
      (21)


      Tabled 1
      ParameterValueUnitSource
      Vf0.782mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.1mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Vr55.83mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,r5mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,AMP0.016mM(
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      )

      Lactate dynamics

      PYR+NADH+HL+NAD
      (R21)


      JLDH=Vf[PYR]Km,f+[PYR]1+Km,NADH[NADH]
      (22)


      Tabled 1
      ParameterValueUnitSource
      Vf23.93mM/s(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,f0.125mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Km,NADH0.001mM(
      • Kashiwaya Y.
      • Sato K.
      • Tsuchiya N.
      • Thomas S.
      • Fell D.A.
      • Veech R.L.
      • Passonneau J.V.
      Control of glucose utilization in working perfused rat heart.
      )
      Le+HeL+H
      (R22)


      JLHX=Vf[L]e[L]e+Km,fVr[L]c[L]c+Km,r
      (23)


      Tabled 1
      ParameterValueUnitSource
      Vf0.048mM/s(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      Km,f2.2mM(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      Vr0.182mM/s(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      Km,r6.92mM(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      L+HLA
      (R23)


      JLA=XLA([L]c[H]cKm[LA]c)
      (24)


      Tabled 1
      ParameterValueUnitSource
      XLA10,0001/s(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )
      Km1.259×107mM(
      • Ch’en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      )

      Other

      Mb+O2MbO2
      (R24)


      JMB=ka[Mb]c[O2]ckd[MbO2]c
      (25)


      Tabled 1
      ParameterValueUnitSource
      ka15,4001/(mMs)(
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      )
      kb601/s(
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      )
      fADP+MgmADP
      (R25)


      βADP=12(KMgADP+[Mg]c,tot+[mADP]cKMgADP+[Mg]c,tot4[mADP]c[Mg]c,tot)
      (26)


      Tabled 1
      ParameterValueUnitSource
      KMgADP0.347mM(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ,
      • Korzeniewski B.
      Theoretical studies on the regulation of oxidative phosphorylation in intact tissues.
      )
      [Mg]c,tot5mM(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )

      Mitochondrial intermembrane space

      fADP+MgmADP
      (R26)


      fATP+MgmATP
      (R27)


      JMgADPi=XJMgA([ADPf]i[Mg]iKMgADP[mADP]i)
      (27)


      JMgATPi=XMgA([fATP]i[Mg]iKMgATP[mATP]i)
      (28)


      Tabled 1
      ParameterValueUnitSource
      XMgA10001/(mMs)(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      )
      KMgADP0.347mM(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ,
      • Korzeniewski B.
      Theoretical studies on the regulation of oxidative phosphorylation in intact tissues.
      )
      KMgATP0.024mM(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.
      ,
      • Korzeniewski B.
      Theoretical studies on the regulation of oxidative phosphorylation in intact tissues.
      )

      Phosphate transporter

      H2PO4HPO4+H
      (R28)


      JPi1=XPi1([H]x[H2PO4]i[H]i[H2PO4]x)[H2PO4]i+kPiHt
      (29)


      Tabled 1
      ParameterValueUnitSource
      XPi13.394 × 1051/s(
      • Beard D.A.
      A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation.