Advertisement

Modeling oxygen requirements in ischemic cardiomyocytes

  • Anthony D. McDougal
    Affiliations
    Departments of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
    Search for articles by this author
  • C.Forbes Dewey Jr.
    Correspondence
    To whom correspondence should be addressed: Dept. of Mechanical Engineering, Massachusetts Institute of Technology, Rm. 3-354, 77 Massachusetts Ave., Cambridge, MA 02139. Tel.: 617-253-2235;
    Affiliations
    Departments of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139

    Biological Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
    Search for articles by this author
Open AccessPublished:May 09, 2017DOI:https://doi.org/10.1074/jbc.M116.751826
      Heart disease remains the leading cause of death globally. Although reperfusion following myocardial ischemia can prevent death by restoring nutrient flow, ischemia/reperfusion injury can cause significant heart damage. The mechanisms that drive ischemia/reperfusion injury are not well understood; currently, few methods can predict the state of the cardiac muscle cell and its metabolic conditions during ischemia. Here, we explored the energetic sustainability of cardiomyocytes, using a model for cellular metabolism to predict the levels of ATP following hypoxia. We modeled glycolytic metabolism with a system of coupled ordinary differential equations describing the individual metabolic reactions within the cardiomyocyte over time. Reduced oxygen levels and ATP consumption rates were simulated to characterize metabolite responses to ischemia. By tracking biochemical species within the cell, our model enables prediction of the cell's condition up to the moment of reperfusion. The simulations revealed a distinct transition between energetically sustainable and unsustainable ATP concentrations for various energetic demands. Our model illustrates how even low oxygen concentrations allow the cell to perform essential functions. We found that the oxygen level required for a sustainable level of ATP increases roughly linearly with the ATP consumption rate. An extracellular O2 concentration of ∼0.007 mm could supply basic energy needs in non-beating cardiomyocytes, suggesting that increased collateral circulation may provide an important source of oxygen to sustain the cardiomyocyte during extended ischemia. Our model provides a time-dependent framework for studying various intervention strategies to change the outcome of reperfusion.

      Introduction

      Globally, ischemic heart disease was the leading cause of death in 2012 and the most rapidly increasing cause of death during 2000–2012 (
      United States Burden of Disease Collaborators
      ,
      World Health Organization
      ). Ischemia occurs when the circulation of the blood is restricted, thereby limiting the delivery of nutrients and removal of metabolic by-products. The heart can be saved by restoring blood flow using reperfusion techniques such as percutaneous coronary intervention (
      • Yellon D.M.
      • Hausenloy D.J.
      Myocardial reperfusion injury.
      ). However, even as the heart is saved, reperfusion carries the risk of damaging additional heart tissue. This damage is termed ischemia/reperfusion injury. Various proposals for novel reperfusion techniques have been advanced, but our understanding of ischemia/reperfusion injury and our ability to mitigate its risk are significantly lacking (
      • Hausenloy D.J.
      • Yellon D.M.
      Myocardial ischemia-reperfusion injury: a neglected therapeutic target.
      • Kalogeris T.
      • Baines C.P.
      • Krenz M.
      • Korthuis R.J.
      Cell biology of ischemia/reperfusion injury.
      ,
      • Schwartz Longacre L.
      • Kloner R.A.
      • Arai A.E.
      • Baines C.P.
      • Bolli R.
      • Braunwald E.
      • Downey J.
      • Gibbons R.J.
      • Gottlieb R.A.
      • Heusch G.
      • Jennings R.B.
      • Lefer D.J.
      • Mentzer R.M.
      • Murphy E.
      • Ovize M.
      • et al.
      New horizons in cardioprotection: recommendations from the 2010 National Heart, Lung, and Blood Institute Workshop.
      ,
      • Ferdinandy P.
      • Hausenloy D.J.
      • Heusch G.
      • Baxter G.F.
      • Schulz R.
      Interaction of risk factors, comorbidities, and comedications with ischemia/reperfusion injury and cardioprotection by preconditioning, postconditioning, and remote conditioning.
      ,
      • Bell R.M.
      • Bøtker H.E.
      • Carr R.D.
      • Davidson S.M.
      • Downey J.M.
      • Dutka D.P.
      • Heusch G.
      • Ibanez B.
      • Macallister R.
      • Stoppe C.
      • Ovize M.
      • Redington A.
      • Walker J.M.
      • Yellon D.M.
      9th Hatter Biannual Meeting: position document on ischaemia/reperfusion injury, conditioning and the ten commandments of cardioprotection.
      • Kalogeris T.
      • Baines C.P.
      • Krenz M.
      • Korthuis R.J.
      Ischemia/reperfusion.
      ).
      It is imperative we identify the quantitative conditions that exist in the cardiomyocyte following a period of ischemia. To date, a number of groups have modeled the dynamics of heart tissue (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      • Cortassa S.
      • Aon M.A.
      • O'Rourke B.
      • Jacques R.
      • Tseng H.-J.
      • Marbán E.
      • Winslow R.L.
      A computational model integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte.
      ,
      • 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.
      ,
      • Korzeniewski B.
      • Noma A.
      • Matsuoka S.
      Regulation of oxidative phosphorylation in intact mammalian heart in vivo.
      ,
      • Tepp K.
      • Timohhina N.
      • Chekulayev V.
      • Shevchuk I.
      • Kaambre T.
      • Saks V.
      Metabolic control analysis of integrated energy metabolism in permeabilized cardiomyocytes–experimental study.
      ,
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      • Yaniv Y.
      • Stanley W.C.
      • Saidel G.M.
      • Cabrera M.E.
      • Landesberg A.
      The role of Ca2+ in coupling cardiac metabolism with regulation of contraction.
      ). However, only Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ) model ischemic cardiomyocytes as fully deprived of oxygen, and their model does not provide quantitative data of the detailed transition in time between a fully oxygenated myocardium and one that is deficient of oxygen. 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.
      ) represent the onset of ischemia from normal conditions (for 30 s). An alternative approach is taken by Karlstädt et al. (
      • Karlstädt A.
      • Fliegner D.
      • Kararigas G.
      • Ruderisch H.S.
      • Regitz-Zagrosek V.
      • Holzhütter H.-G.
      CardioNet: a human metabolic network suited for the study of cardiomyocyte metabolism.
      ) to construct an extensive model that identifies minimum substrate and oxygen requirements for normal function of the cardiomyocyte, but it does not address the changes that develop during hypoxia. The model by Zhou et al. (
      • Zhou L.
      • Salem J.E.
      • Saidel G.M.
      • Stanley W.C.
      • Cabrera M.E.
      Mechanistic model of cardiac energy metabolism predicts localization of glycolysis to cytosolic subdomain during ischemia.
      ) is the most “complete” in regards to including all of the subdomains of metabolism during ischemia, and elements of glycolysis, fatty acid metabolism, the citric acid cycle, and oxidative phosphorylation are included. Zhou et al. (
      • Zhou L.
      • Salem J.E.
      • Saidel G.M.
      • Stanley W.C.
      • Cabrera M.E.
      Mechanistic model of cardiac energy metabolism predicts localization of glycolysis to cytosolic subdomain during ischemia.
      ) use their model to illustrate the transition from oxygenated conditions to partial deoxygenation.
      Here, we build upon these efforts to present a model that describes metabolism through various levels of deoxygenation and metabolic demand. Our aim is to identify the precursor conditions developed during ischemia that leave cells susceptible to ischemia/reperfusion injury. The conditions that are suspected to play a role in ischemia/reperfusion injury include calcium overload, oxidative stress, mitochondrial dysfunction, and signaling cues from death proteins. These conditions develop due to a lack of energy from ATP during the first stage of ischemia/reperfusion injury (ischemia up to the moment of reperfusion) and further evolve during a second stage, which begins upon the first instant of reperfusion. To understand the mechanisms driving ischemia/reperfusion injury, both stages require quantification of two aspects, the energy balance of the cell as well as the presence of biochemical species (and their interactions). We begin a first approach toward these objectives by evaluating the energy production of the cell during ischemia leading up to reperfusion, thereby tracking some of the metabolites as well. More specifically, we ask how a drop in oxygen affects the ATP available for cellular functions. To this end, we have developed a model of cardiomyocyte metabolism that accounts for cytoplasmic metabolism via glycolysis, mitochondrial oxidation of pyruvate, ATP buffering, and ion transport. Glycolysis and glycogenolysis are described in particular detail to account for anaerobic metabolism.
      A significant amount of controversy surrounds the topics of energy demand, production, and regulation in cardiomyocytes (
      • Aon M.A.
      • Cortassa S.
      Mitochondrial network energetics in the heart.
      ). One focal aspect of this debate is the mechanism by which ATP consumption is regulated, whether by calcium, inorganic phosphate, creatine, or otherwise (
      • 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.
      ). Rather than evaluating these possible feedback mechanisms, we focus on understanding what ATP consumption rates are sustainable prior to reperfusion, given the resources and the method of metabolism of the cardiomyocyte. Using this approach, we are able to identify sets of reasonable and unreasonable outcomes for the heart cell by evaluating simulations across the entire range of possible energy demands.
      Our model uses reaction parameters from the literature to simulate the time-development of cellular conditions following the onset of hypoxia. From these representative simulations, we predict how the availability of oxygen impacts the sustainability of a cardiomyocyte's energy demands, and hence survival. This is a crucial step in understanding the changes induced by reperfusion of the tissue, which has clinically been found to induce additional tissue injury.
      In normal heart function, over 95% of ATP regeneration comes from oxidative phosphorylation in the mitochondria; moreover, 50–70% of ATP comes from fatty acids (
      • Lopaschuk G.D.
      • Ussher J.R.
      • Folmes C.D.
      • Jaswal J.S.
      • Stanley W.C.
      Myocardial fatty acid metabolism in health and disease.
      ). However, during ischemia, β-oxidation of fatty acids slows while anaerobic glycolysis increases. One proposed strategy in the face of ischemia and reperfusion is to improve glucose oxidation and its coupling to glycolysis (
      • Fillmore N.
      • Lopaschuk G.D.
      Malonyl CoA: a promising target for the treatment of cardiac disease.
      ). In this vein, we explore the extent to which glycolysis can provide ATP as well as pyruvate for oxidation and leave the contribution and inhibition of fatty acids for another study. Similarly, we do not detail the reactions of the mitochondria, but rather we explore its ensemble performance to find the cell's best-case energy performance under the challenging conditions of ischemia.
      The model contains the following three regions (see Fig. 1): the cytoplasm of the cell, the extracellular space, and a blood vessel. The mitochondria is included within the cell as a single representative overall reaction of oxidative phosphorylation of pyruvate. Each region was treated as a well mixed compartment. Because the purpose of this study was to evaluate the ability of glycolysis to support energy demands, geometry and space considerations were left for future development. Reactions were included to describe glycolysis (as depicted in Fig. 2), glycogenesis/glycogenolysis, mitochondrial oxidation of pyruvate, ATP consumption, ATP buffering, effects of ion transport on pH, and oxygen permeation across compartments. Parameters for all reactions were taken from values reported in experiments, other computational models, or otherwise in the literature. As our model focuses on glycolytic capabilities, we make some assumptions regarding other processes. In real tissue, one would expect changes to NAD, calcium, etc. to influence many of the reactions. In general, we have accounted for the effect of these molecules within the expression of the reaction itself (for example, hexokinase depends on the changing presence ATP). In the case of some molecules (e.g. NAD/NADH), we account for their impact on the reaction velocity but have held the concentration of that species constant due to its interaction with many other processes. Thus our model predominantly describes the development of species within the glycolytic metabolism to first order.
      Figure thumbnail gr1
      Figure 1Overview of the metabolic components included in the model. We focus on the increased role of glycolysis during hypoxia (
      • Lopaschuk G.D.
      • Ussher J.R.
      • Folmes C.D.
      • Jaswal J.S.
      • Stanley W.C.
      Myocardial fatty acid metabolism in health and disease.
      ) and do not explore here the role of fatty acid metabolism. All parameters are taken from the literature. Our model evaluates system reactions using a specified fraction of available oxygen and of the ATP consumption rate.
      Figure thumbnail gr2
      Figure 2Reactions included in the model of glycolysis. Abbreviations for biochemical species are found in . Reactions for glycolysis are shown in orange, and the reactions that the model simulates as interacting with components of glycolysis are shown in black.
      Table 1Species evaluated in the model
      SpeciesAbbreviationInitial valueSource
      mm
      1,3-Bisphosphoglycerate13BPG8.69 × 10−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.
      2-Phosphoglycerate2PG0.00900
      • 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.
      3-Phosphoglycerate3PG0.0710
      • 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.
      Adenosine diphosphateADP0.00
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Adenosine monophosphateAMP1.00 × 10−5
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Initial [AMP] is considered to be 0 by Ch'en et al. (10). Here we have set the initial value to be a very low, negligible quantity to avoid a “divide by zero” error.
      Adenosine triphosphateATP7.00
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      BicarbonateHCO312.5Calculated from Ref.
      • Lagadic-Gossmann D.
      • Buckler K.J.
      • Vaughan-Jones R.D.
      Role of bicarbonate in pH recovery from intracellular acidosis in the guinea-pig ventricular myocyte.
      Bicarbonate, extracellularHCO3,e25.0
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      CreatineCr0.00
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Creatine phosphateCrP25.0
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Dihydroxyacetone phosphateDHAP0.0360
      • 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.
      Fructose 1,6-bisphosphateF16BP6.78 × 10−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.
      Fructose 6-phosphateF6P0.0410
      • 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.
      Glucose 1,6-bisphosphateG16BP0.00700
      • 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.
      Glucose 1-phosphateG1P0.0200
      • 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.
      Glucose 6-phosphateG6P0.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.
      Glyceraldehyde 3-phosphateGAP0.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.
      Glycerol 3-phosphateGly-3-P0.158
      Data were calculated from Denton et al. (66) and Kashiwaya et al. (60).
      GlycogenGlg21.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.
      HydrogenH+7.94 × 10−5
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Inorganic phosphatePi7.00
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      LactateL0.247
      Data were calculated using the value for L−e and equations in Ch'en et al. (10).
      Lactic acidLaH156
      Data were calculated using the value for L−e and equations in Ch'en et al. (10).
      Myoglobin, freeMb0.00543
      Endeward (65) identifies [Mbtot] = 0.19 mm. The initial values for [Mb] and [MbO] are the resulting equilibrium values given the association rates for oxygen binding (see below) and [O2].
      Myoglobin, saturatedMbO20.185
      Endeward (65) identifies [Mbtot] = 0.19 mm. The initial values for [Mb] and [MbO] are the resulting equilibrium values given the association rates for oxygen binding (see below) and [O2].
      Oxygen (intracellular)O20.110
      Early iterations of the model were run with the given [O2, v] value, and observed to find the steady [O2, e] and [O2].
      Oxygen, extracellularO2,e0.129
      Early iterations of the model were run with the given [O2, v] value, and observed to find the steady [O2, e] and [O2].
      PhosphoenolpyruvatePEP0.0130
      • 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.
      PyruvatePYR0.0550
      • 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.
      Uridine diphosphateUDP0.0340
      • 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.
      Uridine diphosphate glucoseUDPG0.0990
      • 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.
      Uridine triphosphateUTP0.194
      • 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.
      a Initial [AMP] is considered to be 0 by Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ). Here we have set the initial value to be a very low, negligible quantity to avoid a “divide by zero” error.
      b Data were calculated from Denton et al. (
      • Denton R.M.
      • Yorke R.E.
      • Randle P.J.
      Measurement of concentrations of metabolites in adipose tissue and effects of insulin, alloxan-diabetes and adrenaline.
      ) and 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.
      ).
      c Data were calculated using the value for Le and equations in Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ).
      d Endeward (
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      ) identifies [Mbtot] = 0.19 mm. The initial values for [Mb] and [MbO] are the resulting equilibrium values given the association rates for oxygen binding (see below) and [O2].
      e Early iterations of the model were run with the given [O2, v] value, and observed to find the steady [O2, e] and [O2].
      Consistent with our goal to find how the availability of oxygen impacts the sustainability of particular energy demands in the myocyte, our model uses a fixed ATP consumption rate rather than incorporating one of many debated mechanisms of feedback regulation. Instead of having a deterministic ATP consumption, here we map out the cellular response by changing only two parameters from simulation to simulation: the oxygen level present in the ischemic blood vessel and the level of ATP consumption. In every simulation, the initial blood vessel oxygen was set to 0.133 mm (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ) and subsequently dropped to a fixed ischemic level after 100 s. Simulations continued for 5000 s, an interval that allowed all cases to reestablish a steady state. Over 200 trials were simulated in CellDesigner with time steps of 5 s. The details of the model are presented under “Computational methods.”

      Results

      The simulations predict the metabolic behavior of the cardiomyocyte upon hypoxia. Fig. 3 represents typical time histories of ATP, ADP, and AMP and the energy buffers creatine phosphate and glycogen following a drop in blood oxygen. The concentration of ATP is maintained by a sequence of three energy buffers: we observe a drop in creatine phosphate, followed by the fall of glycogen, and finally a rise in AMP as ADP is sacrificed. Throughout these events, the reactions of glycolysis are providing an anaerobic source of ATP. Fig. 3A illustrates a typical response to deoxygenation after a brief period of stable, oxygenated conditions. After 100 s of simulation, the supplied oxygen is decreased and held at a new low level. After the energy buffers are depleted, [ATP] falls to zero in about 550 s, at which point the consumption rate of ATP is no longer sustainable and ATP hydrolysis fails. In comparison, Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ) used a completely anaerobic model and found [ATP] to last roughly 800 s. 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.
      ) provide another point of comparison regarding creatine phosphate (CrP)
      The abbreviations used are: CrP, creatine phosphate; ROS, reactive oxygen species; SOD, superoxide dismutase; SBML, Systems Biology Markup Language; G6P, Glc-6-P; abbreviations used for biochemical species are found in Table 1.
      : when totally occluding canine coronary blood flow, CrP fell to 10 mm in 35 s; our simulations take 110 s to show the same change but under incomplete occlusion (the simulations of 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.
      ) are not extended to full depletion of CrP).
      Figure thumbnail gr3
      Figure 3Metabolite response to deoxygenation. Two example simulations of metabolite response to a decrease in the supply of oxygen in the blood. A, at the outset, oxygen in the blood [O2, V] = 0.133 mm and the system stays at steady state; after 100 s (arrow), [O2, V] falls to 0.0150 mm. B, [O2, V] falls from 0.133 to 0.0175 mm after 100 s (arrow). For both simulations, rATP = 100%, i.e. constant ATP consumption rate of 0.570 mm/s, indicative of a beating heart in a resting subject. Shown are predicted concentrations of: creatine phosphate (light blue), glycogen (dark blue), ATP (orange), ADP (yellow), and AMP (green). The arrow indicates the decrease of oxygen in the blood at 100 s.
      It should be reiterated that our simulations drive ATP consumption at a specified rate, rather than being controlled by a particular feedback mechanism. In the real cardiomyocyte, as resources are exhausted and by-products build up, this feedback will work to slow the consumption of ATP, prolonging the life of the cell. Thus, Fig. 3A can serve as a worst case timeline for the depletion of energy buffers.
      If the oxygen supply is decreased to an even lower value, the rate of mitochondrial oxidation of pyruvate is additionally decreased, further hampering the greatest energy source. Consequently, the energy buffers will deplete faster in an effort to keep [ATP] at high levels, and the cell will crash more rapidly. Conversely, if the oxygen is slightly increased, mitochondrial reactions are able to provide more energy to supplement anaerobic metabolism, which lessens the pull on the buffers and extends the survival of the cell (see Fig. 3B). (A variable ATP consumption rate will alter the time to ATP depletion; this is explored under “Discussion.”)
      After the transient behavior of the system, we see that the final outcome of each simulation is determined by the amount of oxygen available and the amount of ATP being consumed. If the supplied oxygen level is sufficiently high, the concentration of ATP shows no perceptible change (see Fig. 4A). However, there is a transition range of oxygen concentrations that lead to an intermediate steady-state level of [ATP] (as in Fig. 3B). This result is surprising because it suggests that although ATP demand can be matched by ATP regeneration, the cell does not do so until after consuming its backup supply and allowing some decrease in [ATP]. Thus we see that, for a specified ATP consumption rate, lowering oxygen can affect the steady levels of ATP, resulting in catastrophically low concentrations at O2 levels below a fairly sharp threshold. As one may expect, if we decrease the energy consumption rate under low O2 conditions, [ATP] no longer plummets (Fig. 4B).
      Figure thumbnail gr4
      Figure 4[ATP] trends for various simulation parameters. A, comparison of [ATP] trends for constant ATP consumption for a beating cell (0.570 mm/s) upon decreasing blood vessel oxygen from initial [O2, V] = 0.133 mm. Solid line, [O2, V] falls to 0.0150 mm, identical to A. Dashed line, [O2, V] falls to 0.0175 mm, identical to B. Dotted line, [O2, V] falls to 0.0200 mm. The arrow indicates the decrease of oxygen in the blood at 100 s for all cases. B, comparison of [ATP] trends when blood vessel oxygen falls from initial [O2, V] = 0.133 to 0.0150 mm, for various ATP consumption rates. Solid line, rATP = 100%, identical to A and A. Dot-dashed line, rATP = 90%. Dashed line, rATP = 82%. Dotted line, rATP = 80%. The arrow indicates the decrease of oxygen in the blood at 100 s for all cases.
      This leads to the following question: when is [ATP] too low to sustain cell energy demands? Even if there persists an extremely small, steady amount of [ATP] for energy turnover, any further perturbation could lead to decreasing [ATP] such that the ATP consumption rate cannot be sustained. For the purposes of our computations, we treat 0.1% of the initial [ATP] (= 0.007 mm) as the lowest sustainable concentration of ATP. Similarly, we define “imperceptible changes” of [ATP] as less than a 0.1% change of the original concentration of 7.000 mm. At the lower bound, we recognize that concentrations may be too low to model using reaction rates, and the biochemical interactions are better represented by the Langevin equation and become stochastic in nature. At such low levels, the molecular interactions will be slower than the reaction rates predict.
      We iterated our model to find all the combinations of [O2] and ATP consumption rates that lead to these [ATP] thresholds. The ATP consumption rates from these simulations and the resulting extracellular [O2, e] are plotted in Fig. 5.
      Figure thumbnail gr5
      Figure 5Curves of steady-state ATP concentration, arising from ATP consumption and extracellular oxygen. Simulations were run over all ATP consumption rates and oxygen concentrations. Although the blood vessel oxygen continues to be the simulation input as in , for this figure we plot the extracellular oxygen to emphasize the oxygen available in the tissue. The drop in oxygen from the vessel to the extracellular space will vary depending on the ATP consumption rate. Solid red line, those consumption rates and [O2, V] that lead to a steady state of [ATP] equal to 99.9% of the original ATP. Solid blue line, those consumption rates and [O2, V] that lead to a steady state of [ATP] equal to 0.1% of the original ATP. Dashed red and blue lines, range of conditions that lead to intermediate [ATP] at steady state.
      In the right-hand region of Fig. 5, [ATP] shows practically no change from the original levels, and so from a metabolic standpoint, the oxygen supply and energy demand result in a sustainable state for the tissue. In the left-hand region, [ATP] is at such low levels that even the slightest increase of ATP consumption or the slightest decrease of O2 will cause [ATP] to be insufficient. In the middle region, we see a thin range of oxygen concentrations that represent the intermediate transition from sustainable to unsustainable conditions. This “metabolic twilight” region neither exhausts ATP nor maintains it at normal levels.
      These three distinct regions of metabolic outcomes offer a new and useful framework for interpreting cell performance under hypoxic conditions. Yet all of these changes occur only once the oxygen supply is low, under 0.0160 mm extracellular O2 (about 12.8 mm Hg). This falls just under the venous PO2 values of 15–20 mm Hg reported in the literature (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ,
      • Rumsey W.L.
      • Pawlowski M.
      • Lejavardi N.
      • Wilson D.F.
      Oxygen pressure distribution in the heart in vivo and evaluation of the ischemic “border zone.”.
      ). Additionally, the fact that such small changes in oxygen can result in drastically different metabolic outcomes is consistent with the well known fact that patients with significant collateral circulation upon X-ray angiography have improved outcomes following reperfusion (
      • Seiler C.
      • Stoller M.
      • Pitt B.
      • Meier P.
      The human coronary collateral circulation: development and clinical importance.
      ).
      Given this map, we can begin to identify some parameters that are key for the survival of the cardiomyocyte. Das and Harris (
      • Das A.M.
      • Harris D.A.
      Control of mitochondrial ATP synthase in heart cells: inactive to active transitions caused by beating or positive inotropic agents.
      ) analyzed the energy consumed by rat cardiomyocytes under stimulation. Unstimulated, non-contracting cells were found to consume 3.7 μmol ATP/min/mg. Cells that were stimulated were found to have an increasing amount of ATP consumption, plateauing at 6.6 μmol ATP/min/mg. Thus, Das and Harris (
      • Das A.M.
      • Harris D.A.
      Control of mitochondrial ATP synthase in heart cells: inactive to active transitions caused by beating or positive inotropic agents.
      ) indicate that a non-contracting cell consumes 56% of the energy that a contracting cell consumes (see Fig. 6).
      Figure thumbnail gr6
      Figure 6Changes to energy consumed by the rat heart. Summary of data found in Das and Harris (
      • Das A.M.
      • Harris D.A.
      Control of mitochondrial ATP synthase in heart cells: inactive to active transitions caused by beating or positive inotropic agents.
      ). Non-contracting oxygenated cells were found to expend an additional 79% of energy when stimulated to contract.
      Rolfe and Brown (
      • Rolfe D.F.
      • Brown G.C.
      Cellular energy utilization and molecular origin of standard metabolic rate in mammals.
      ) provide a review investigating the allocation of energy among various cellular processes. By compiling their data as listed in Fig. 7, we can consider that 79% of energy is directed toward contraction, whereas 21% is directed toward “cellular maintenance.” From these estimates, we can define the range of oxygen that might possibly sustain cardiomyocytes, allowing us to determine the cells that are targeted during reperfusion.
      Figure thumbnail gr7
      Figure 7Percentage of energy consumed by processes in the rat heart. Summary of data presented in Rolfe and Brown (
      • Rolfe D.F.
      • Brown G.C.
      Cellular energy utilization and molecular origin of standard metabolic rate in mammals.
      ). The values reported by Rolfe and Brown are given as ranges; the values shown here are a proportional distribution of those values. They estimate that activities associated with contraction make up about roughly 79% of energy consumption, whereas the activities necessary for cellular maintenance lead to 21%.
      Combining the ATP state information in Fig. 5 with the ATP consumption data in Figure 6, Figure 7 to create Fig. 8, we find that when the extracellular O2 drops below 0.011 mm extracellular oxygen in a non-beating heart, ATP begins to fall from its sustainable state (red line in Fig. 8). In a heart consuming energy only for maintenance purposes (21%), the oxygen requirement for unaltered ATP (Fig. 8, red line) is around 0.006 mm. For maintenance purposes, the unsustainable region (Fig. 8, blue line) is reached at about 0.004 mm extracellular oxygen.
      Figure thumbnail gr8
      Figure 8Steady-state ATP concentration curves compared with cellular energy requirements. A comparison of the results from Figure 5, Figure 6, Figure 7. Top horizontal line, energy consumed by non-beating cardiomyocytes as a percentage of energy required to sustain a beating heart (
      • Das A.M.
      • Harris D.A.
      Control of mitochondrial ATP synthase in heart cells: inactive to active transitions caused by beating or positive inotropic agents.
      ). Bottom horizontal line, energy consumed by key maintenance functions in the cell, as a percentage of energy required to sustain a beating heart (
      • Rolfe D.F.
      • Brown G.C.
      Cellular energy utilization and molecular origin of standard metabolic rate in mammals.
      ). By comparing the results from the model simulations with the energy consumed by non-beating cardiomyocyte functions, we can estimate that a non-beating heart is sustained by extracellular oxygen on the order of 0.005–0.010 mm. Dots, sample results for decreased glucose (0.191 mm): two black dots (left) indicate steady [ATP] of less than 0.007 mm; two green dots (right) indicate steady [ATP] greater than 6.993 mm; and two yellow dots (middle) indicate steady [ATP] between 0.007 and 6.993 mm.
      Of course, ischemia will also produce a decrease in glucose, which will impact the cell's metabolism. Our model holds the glucose level constant during any individual case. The full characterization of the interplay between glucose and oxygen is left for future study. Here, we illustrate some sample outcomes of our model when the amount of available intracellular glucose is moved from 1.91 to 0.191 mm, shown as dots in Fig. 8. For these values of glucose, at ATP consumption rates seen in the non-beating heart, we find that the transition region is nearly unchanged.
      Overall, these simulation results suggest that the cell can balance energy demands using different reactions to different extents. With less oxygen, mitochondrial oxidation of pyruvate is slowed, decreasing the availability of ATP. Over the next 15 min or so, the buffers all attempt to keep the adenosine pool shifted toward ATP as much as possible. Interestingly, the cell uses its buffers despite the fact that, under some conditions, it is able to match ATP supply and demand after the buffers have been depleted. In effect, it appears that the cell's metabolism has multiple safeguards that can be called upon, not only to keep [ATP] at its normal level but also to match ATP supply and demand despite low [ATP].

      Discussion

      The human heart and the cardiomyocyte are incredibly resilient. If an artery is blocked and insufficient O2 is available to keep the ATP at a sustainable level, the myocytes must first drop the energetic load of contraction and then continue to shut the processes down until one reaches a bare minimum that is consistent with the ATP production required to remain viable. If that is not possible, the ATP level rapidly falls to essentially zero, and there is no available energy to take care of critical cellular processes.
      One conclusion of our model is that the cardiac myocyte is able to sustain ATP at a substantial and physiologically useful level for some time following the imposition of a low extracellular O2 level. To achieve this, the cell undergoes a transition where initially the cell's stores of creatine phosphate and glycogen are sacrificed to allow the ATP level to be maintained. The level of energy reserves will persist for times on the order of 90 min if the residual supply of oxygen is sufficiently large. If the residual O2 is eliminated, the energy storage is rapidly consumed, and the ATP and ADP are exhausted in a matter of 10 min. For comparison, the recommended time for “first-medical-contact to device” intervention for patients with ST-elevation myocardial infarction is under 90 min (
      • O'Gara P.T.
      • Kushner F.G.
      • Ascheim D.D.
      • Casey Jr., D.E.
      • Chung M.K.
      • de Lemos J.A.
      • de Ettinger S.M.
      • Fang J.C.
      • Fesmire F.M.
      • Franklin B.A.
      • Granger C.B.
      • Krumholz H.M.
      • Linderbaum J.A.
      • Morrow D.A.
      • Newby L.K.
      • et al.
      2013 ACCF/AHA guideline for the management of ST-elevation myocardial infarction: executive summary.
      ).
      The difference between residual oxygen levels that support physiological concentrations of ATP and residual levels that fail to provide physiological levels of [ATP] is surprisingly small. Depending on the ATP consumption rate of the cell, this difference could be as little as 10% of the “critical level” of O2 where no change in [ATP] is detected. This means that ischemic tissue with a gradient in O2 driven by diffusion from a distant source could produce a spatial change in the tissue, with the tissue nearest the collateral source being salvageable by reperfusion and the tissue slightly further away being at risk. This is schematically illustrated in Fig. 9. The extent of this “twilight range” depends on the external [O2] gradient.
      Figure thumbnail gr9
      Figure 9Hypothetical diffusion-reaction curve showing the one-dimensional decay of oxygen as a function of distance from an oxygen source. At lower levels, oxygen decays less rapidly. Consequently, even over significant distances in the tissue we may see small variations in the extracellular oxygen concentration.
      Others (
      • Murphy E.
      • Steenbergen C.
      Mechanisms underlying acute protection from cardiac ischemia-reperfusion injury.
      ,
      • Neely J.R.
      • Morgan H.E.
      Relationship between carbohydrate and lipid metabolism and the energy balance of heart muscle.
      ) have suggested that low ATP is not correlated with cell death. However, the conditions for these claims present [ATP] still within the intermediate “twilight” range, rather than a truly “unsustainable” concentration. For example, early studies from Neely and Morgan (
      • Neely J.R.
      • Morgan H.E.
      Relationship between carbohydrate and lipid metabolism and the energy balance of heart muscle.
      ) observed [ATP] as low as around 3% of normal, compared with our value of 0.1%. At the lower value, slightly increased demand for energy becomes too much for ATP regeneration to keep pace, resulting in a breakdown of the cell's necessary functions (
      • Murphy E.
      • Steenbergen C.
      Mechanisms underlying acute protection from cardiac ischemia-reperfusion injury.
      ).
      The spatial extent of cells transitioning from high to low [ATP] can occur over a range of 0.0015 mm O2. Thus, even slight variations in oxygen that may arise due to extracellular and intracellular oxygen fluctuations may lead to substantial differences in the metabolic viability of a cell. We suggest that it is primarily the cells in this “transition region” that are particularly susceptible to perturbations such as reperfusion. At first glance, the range of [O2] that can effect an intermediate level of ATP seems to imply a remarkably thin tissue layer that is at risk. Many studies have suggested that, indeed, the spatial transition from ischemic to functional tissue is sharp; yet historically, gathering data on transient intermediate states has been difficult (
      • Hearse D.J.
      • Yellon D.M.
      The “border zone” in evolving myocardial infarction: controversy or confusion?.
      ). To illustrate how oxygen may vary spatially, locations along the reaction-diffusion curve have an increasingly smaller gradient of oxygen the further they are from the source, after an initial region of sharp decay (see Fig. 9 for an illustrative drawing). Geometrical considerations will impact oxygen availability across tissue, affecting different regions of the heart, as well as within the cell and across the mitochondrial membranes (
      • Jones D.P.
      • Mason H.S.
      Gradients of O2 concentration in hepatocytes.
      ).
      During ischemia in vivo, the source of oxygen from neighboring tissue may be augmented by anastomosis. This will shift the diffusion curve in Fig. 9 to the right, increasing the amount of tissue with a favorable prognosis. Humans have among the highest attainable coronary collateral flow among mammals, which is key in clinical cases (
      • Seiler C.
      • Stoller M.
      • Pitt B.
      • Meier P.
      The human coronary collateral circulation: development and clinical importance.
      ). Our simulations serve to further emphasize the importance of collateral circulation by demonstrating how even miniscule increases in the amount of oxygen can improve cell energy sustainability (as indicated by higher ATP production). If the oxygen being supplied by collateral circulation is sufficiently high, cardiomyocytes may even be able to sustain a high [ATP] threshold and thereby reduce their risk of ischemia/reperfusion injury. It may be worthwhile to point out that animals generally have lower collateral flow, and thus they may have tissue skewed toward a higher risk of ischemia/reperfusion injury compared with human tissue.
      The computational strategy we adopted was to fix the ATP consumption rate for each simulation case, rather than employ a variable ATP consumption rate as would be expected in real ischemia. Thus, the computed time to ATP level stabilization in our model may differ compared with the corresponding physiological case. If one were to represent ATP consumption as a function of metabolites in the cytoplasm, as is commonly done, one could deterministically simulate the cell's final state for a given amount of oxygen. This outcome will necessarily vary according to the model chosen for ATP consumption. However, we fixed the ATP hydrolysis rate to a different level in each simulation, enabling us to map all possibilities of how both the energy demand and the oxygen availability impact the final steady concentrations of the cell. These results are given in Figure 5, Figure 8.
      Although this approach produces some uncertainty in the time to ATP depletion during the initial transient following ischemia, the final steady state of the cardiomyocyte should be accurately determined by the two variables of ATP consumption rate and local oxygen concentration. An exception to this would be events with major excursions of cellular energy demand, e.g. transient contraction demands on compromised myocytes. Popel (
      • Popel A.S.
      Theory of oxygen transport to tissue.
      ) and others have described how to solve the reaction-diffusion equation for various spatial contexts to track cell metabolism. Our simulations along with the Popel model could be employed as the basis for numerical integration of the oxygen consumption in time and space in a compromised vasculature.
      The lack of O2 and the shortage of ATP will lead to reactions that are involved in cellular death or cellular protection. One factor that is key in considering ischemia/reperfusion injury is the rapid generation of reactive oxygen species (ROS) (
      • Kalogeris T.
      • Baines C.P.
      • Krenz M.
      • Korthuis R.J.
      Cell biology of ischemia/reperfusion injury.
      ,
      • Murphy E.
      • Steenbergen C.
      Mechanisms underlying acute protection from cardiac ischemia-reperfusion injury.
      ,
      • Kim J.-S.
      • Jin Y.
      • Lemasters J.J.
      Reactive oxygen species, but not Ca2+ overloading, trigger pH- and mitochondrial permeability transition-dependent death of adult rat myocytes after ischemia-reperfusion.
      ). Recent studies (
      • Chouchani E.T.
      • Methner C.
      • Nadtochiy S.M.
      • Logan A.
      • Pell V.R.
      • Ding S.
      • James A.M.
      • Cochemé H.M.
      • Reinhold J.
      • Lilley K.S.
      • Partridge L.
      • Fearnley I.M.
      • Robinson A.J.
      • Hartley R.C.
      • Smith R.A.
      • Krieg T.
      • Brookes P.S.
      • Murphy M.P.
      Cardioprotection by S-nitrosation of a cysteine switch on mitochondrial complex I.
      ) focus on the condition of the electron transport chain as the cause for increased ROS. Succinate buildup in response to insufficient ATP production is reported to cause reversal of electron transport and ROS production (
      • Chouchani E.T.
      • Pell V.R.
      • Gaude E.
      • Aksentijević D.
      • Sundier S.Y.
      • Robb E.L.
      • Logan A.
      • Nadtochiy S.M.
      • Ord E.N.
      • 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.
      ). From our simulations, we can predict the cell ATP levels for a given oxygen availability and ATP consumption rate just prior to reperfusion. There is also experimental evidence (
      • Zhang P.
      • Lu Y.
      • Yu D.
      • Zhang D.
      • Hu W.
      TRAP1 provides protection against myocardial ischemia-reperfusion injury by ameliorating mitochondrial dysfunction.
      ) that there is a very substantial (45%) decrease of ATP following reperfusion.
      Weakened cardiomyocytes may have less resources to handle ROS loads (
      • Ferrari R.
      • Guardigli G.
      • Mele D.
      • Percoco G.F.
      • Ceconi C.
      • Curello S.
      Oxidative stress during myocardial ischaemia and heart failure.
      ). The role of antioxidants, such as superoxide dismutase (SOD), during reperfusion remains unclear. Heart outcomes appear to improve with increased activity of SOD (
      • French J.P.
      • Quindry J.C.
      • Falk D.J.
      • Staib J.L.
      • Lee Y.
      • Wang K.K.
      • Powers S.K.
      Ischemia-reperfusion-induced calpain activation and SERCA2a degradation are attenuated by exercise training and calpain inhibition.
      ,
      • Suzuki K.
      • Murtuza B.
      • Sammut I.A.
      • Latif N.
      • Jayakumar J.
      • Smolenski R.T.
      • Kaneda Y.
      • Sawa Y.
      • Matsuda H.
      • Yacoub M.H.
      Heat shock protein 72 enhances manganese superoxide dismutase activity during myocardial ischemia-reperfusion injury, associated with mitochondrial protection and apoptosis reduction.
      ) as well as with administered antioxidants (
      • Schriewer J.M.
      • Peek C.B.
      • Bass J.
      • Schumacker P.T.
      ROS-mediated PARP activity undermines mitochondrial function after permeability transition pore opening during myocardial ischemia–reperfusion.
      ); yet SOD appears to degrade quickly when administered intravenously (
      • Chan P.H.
      Role of oxidants in ischemic brain damage.
      ), and the level of SOD activity throughout reperfusion remains unclear (
      • Dhalla N.S.
      • Elmoselhi A.B.
      • Hata T.
      • Makino N.
      Status of myocardial antioxidants in ischemia–reperfusion injury.
      ). In cells with a very low supply of ATP upon reperfusion, it seems feasible that the absence of ATP available for enzyme production may diminish the activity of SOD, allowing ROS formation to advance unchecked (
      • Fridovich I.
      ).
      Another important player in both ischemia and ischemia/reperfusion injury is the calcium gradient across the inner membrane of the mitochondria. For example, Ca2+-ATPase pumps will weaken as ATP is decreased (
      • Glancy B.
      • Balaban R.S.
      Role of mitochondrial Ca2+ in the regulation of cellular energetics.
      ). In turn, the calcium reduces the amount of ATP that can be produced, effectively dismantling the cell's energy production and accelerating the death of the cell. The rise of calcium may also trigger the opening of the mitochondrial permeability transition pore (
      • Di Lisa F.
      • Canton M.
      • Carpi A.
      • Kaludercic N.
      • Menabò R.
      • Menazza S.
      • Semenzato M.
      Mitochondrial injury and protection in ischemic pre- and postconditioning.
      ). Our model shows that a small increase in the amount of oxygen available from collateral circulation may dramatically improve the availability of ATP, which could improve the handling of Ca2+ and ultimately improve the outcome of the cell.
      The model presented here would benefit from additional development. As mentioned previously, a major challenge in modeling cardiomyocytes arises from the debate on how to specify ATP consumption rates as functions of metabolite concentrations (
      • 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.
      ,
      • Karlstädt A.
      • Fliegner D.
      • Kararigas G.
      • Ruderisch H.S.
      • Regitz-Zagrosek V.
      • Holzhütter H.-G.
      CardioNet: a human metabolic network suited for the study of cardiomyocyte metabolism.
      ). Unfortunately, there is very little agreement to guide this modeling.
      We have set the effect of some molecular concentrations to be constant during the onset of ischemia, even though they are certainly changing. Examples include NADH, pH, glucose, and intermediate products of glycolysis that participate in other reactions (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ,
      • 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.
      ,
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      ,
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      ,
      • Vinnakota K.
      • Kemp M.L.
      • Kushmerick M.J.
      Dynamics of muscle glycogenolysis modeled with pH time course computation and pH-dependent reaction equilibria and enzyme kinetics.
      ). One illustration of the effects due to these changing concentrations is demonstrated by Vinnakota et al. (
      • Vinnakota K.
      • Kemp M.L.
      • Kushmerick M.J.
      Dynamics of muscle glycogenolysis modeled with pH time course computation and pH-dependent reaction equilibria and enzyme kinetics.
      ) as they explore the role of pH in the context of glycogenolysis. They show how accounting for variation of the pH does affect the final buildup of by-products by around 30%, while producing similar intermediate results.
      Our model predicts that changing the extracellular glucose concentration by a factor of 10 does not substantially change the history of ATP levels in the cell immediately following the onset of ischemia (see Fig. 8). In the same vein, a sensitivity analysis of various parameters in the model may give further insight regarding the robustness of various functions in the cell. Sensitivity analysis could also indicate the degree to which parameters, such as Km, O2 for oxidative phosphorylation, are sufficient or warrant further refinement.
      Additionally, although our model offers detail on glycolysis, other models emphasize details in other aspects of metabolism, such as the tricarboxylic acid cycle and the electron transport chain. Further details of all the relevant reactions within the mitochondria remain an active area of research, as summarized by a recent review paper (
      • Murphy E.
      • Ardehali H.
      • Balaban R.S.
      • DiLisa F.
      • Dorn 2nd., G.W.
      • Kitsis R.N.
      • Otsu K.
      • Ping P.
      • Rizzuto R.
      • Sack M.N.
      • Wallace D.
      • Youle R.J.
      American Heart Association Council on Basic Cardiovascular Sciences, Council on Clinical CardiologyCouncil on Functional Genomics and Translational Biology
      Mitochondrial function, biology, and role in disease.
      ). These include the following: elements of mitochondrial structure, such as the membrane permeability transition pore (
      • Bell C.J.
      • Bright N.A.
      • Rutter G.A.
      • Griffiths E.J.
      ATP regulation in adult rat cardiomyocytes time-resolved decoding of rapid mitochondrial calcium spiking imaged with targeted photoproteins.
      ); ion regulation and transport, including mitochondrial potassium and calcium (
      • Schulz R.
      • Di Lisa F.
      Mitochondrial potassium homeostasis: a central player in cardioprotection.
      ); pools of key molecular species, for example NAD/NADH (
      • Moxley M.A.
      • Beard D.A.
      • Bazil J.N.
      Global kinetic analysis of mammalian E3 reveals pH-dependent NAD+/NADH regulation, physiological kinetic reversibility, and catalytic optimum.
      ); and information on molecular behavior, such as the details regarding ROS generation (
      • 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.
      ,
      • Gauthier L.D.
      • Greenstein J.L.
      • O'Rourke B.
      • Winslow R.L.
      An integrated mitochondrial ROS production and scavenging model: implications for heart failure.
      ). A complete understanding of metabolism during extreme conditions will depend on simultaneous integration of the strengths of each of these models. The integration process and the additional computational complexity can be carried out with reasonable effort (
      • Ayyadurai V.A.
      • Dewey C.F.
      CytoSolve: a scalable computational method for dynamic integration of multiple molecular pathway models.
      ). For example, the mitochondrial model iAS253 for flux balance analysis by Smith and Robinson (
      • Smith A.C.
      • Robinson A.J.
      A metabolic model of the mitochondrion and its use in modelling diseases of the tricarboxylic acid cycle.
      ) was employed by 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.
      • 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.
      ) to demonstrate how the buildup of succinate during ischemia may drive reverse electron transport and the generation of ROS during reperfusion. Future integration of the model iAS253, our model, and spatial determination of oxygen in cardiac tissue may enable us to predict the entire time course of ischemia/reperfusion injury.
      Looking ahead to extending our current model into predictions following reperfusion, we are very interested in seeing what happens not only to the ROS but also in what happens to anaerobic and aerobic ATP production. The transient state of the cardiac myocyte is key. One cannot understand this process using equilibrium conditions where nothing is changing because the transient in ATP following reperfusion is a major factor. If ATP is decreased significantly, even if only for a short time, and if the influx of metabolites and oxygen during reperfusion results in non-optimal ATP consumption, we could have a serious energy crisis that would put the myocyte in a very precarious position prior to the possibility of replenishing the ATP levels through large-scale ATP production in the mitochondria. The initial rush of reperfusion may transiently drop the ATP/ADP levels sufficiently low that irreversible damage occurs to the myocyte.

      Computational methods

      We created a computational model of a cardiomyocyte using the Systems Biology Markup Language (SBML) (
      • Finney A.
      • Hucka M.
      Systems biology markup language: Level 2 and beyond.
      ,
      • Hucka M.
      • Finney A.
      • Sauro H.M.
      • Bolouri H.
      • Doyle J.C.
      • Kitano H.
      • Arkin A.P.
      • Bornstein B.J.
      • Bray D.
      • Cornish-Bowden A.
      • Cuellar A.A.
      • Dronov S.
      • Gilles E.D.
      • Ginkel M.
      • Gor V.
      • et al.
      The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models.
      ). Reactions were represented as a system of coupled ordinary differential equations and then evaluated in the CellDesigner Simulation environment (
      • Funahashi A.
      • Matsuoka Y.
      • Jouraku A.
      • Morohashi M.
      • Kikuchi N.
      • Kitano H.
      CellDesigner 3.5: a versatile modeling tool for biochemical networks.
      ).

      Biochemical species

      The model includes 56 reactions to explicitly calculate 32 biochemical species while including the effect of other species as detailed in Table 1, Table 2, Table 3. The concentration of every species was recorded for all time steps. Unless otherwise indicated, concentrations refer to values in the cytoplasm. Extracellular concentrations are denoted with a subscript “e” and concentrations in the blood vessel are denoted with a subscript “v.” We record all parameters with an accuracy of three significant figures, regardless of whether these reaction rates were originally reported with more or less significant figures.
      Table 2pH parameters used in the model
      AttributeSymbolSource
      pH, cellularpH 7.10 (initial, variable)
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      pH, extracellularpHe 7.40 (fixed)
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Table 3Species of fixed value in the model
      Fixed speciesSpeciesValueSource
      mm
      Oxidized nicotinamide adenine dinucleotideNAD2.26
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      ,
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      The following data are from Beard (24) and later Wu et al. (15): [NADtot] = 2.97 = [NAD] + [NADH]. Simulations found in Wu et al. (15) indicate that [NADH]/[NAD] ∼0.3 as the heart is actively regenerating ATP.
      Nicotinamide adenine dinucleotideNADH0.714
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      ,
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      The following data are from Beard (24) and later Wu et al. (15): [NADtot] = 2.97 = [NAD] + [NADH]. Simulations found in Wu et al. (15) indicate that [NADH]/[NAD] ∼0.3 as the heart is actively regenerating ATP.
      Lactate, extracellularLe0.330
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      Glucose, intracellular
      We hold the concentration of intracellular glucose constant in this study to focus on the effects of hypoxia. In an ischemic event, glucose could also become depleted over time. We use the concentration of intracellular glucose that was reported as the result of plentiful glucose (10 mm) delivered to the perfuse rat heart (60).
      Glc1.91
      • 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.
      a The following data are from Beard (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ) and later 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.
      ): [NADtot] = 2.97 = [NAD] + [NADH]. Simulations found in 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.
      ) indicate that [NADH]/[NAD] ∼0.3 as the heart is actively regenerating ATP.
      b We hold the concentration of intracellular glucose constant in this study to focus on the effects of hypoxia. In an ischemic event, glucose could also become depleted over time. We use the concentration of intracellular glucose that was reported as the result of plentiful glucose (10 mm) delivered to the perfuse rat heart (
      • 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.
      ).

      Glycolysis

      Fig. 2 depicts the biochemical species and reactions involved in glycolysis. All reaction values of the glycolysis pathway are taken from rat cardiomyocyte data 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.
      ). In our model, reactions are generally represented as reversible Michaelis-Menten equations with modifiers from secondary reactants and products.
      Hexokinase (HK) (Table 4, Reaction 1, and Equation 1):
      Glc+ATPG6P+ADP+H+Reaction 1


      vHK=Vf,HK*Glc*ATPKmf,HKATP+KHKm,ATP-Vr,HK*G6PKmr,HK1+GlcKmf,hk+G6PKmr,HK
      (Eq. 1)


      Phosphoglucose isomerase (PGI) (Table 5, Reaction 2, and Equation 2):
      G6PF6PReaction 2


      vPGI=Vf,PGI*G6PKmf,PGI-Vr,PGI*F6PKmr,PGI1+G6PKmf,PGI+F6PKmr,PGI
      (Eq. 2)


      Phosphofructokinase (PFK) (Table 6, Reaction 3, and Equation 3):
      F6P+ATPF16BP+ADP+H+Reaction 3


      vPFK=Vf,PFK*F6PKmf,PFK+F6P1+KPFKm,ATPATP
      (Eq. 3)


      Fructose-bisphosphate aldolase (F16BP) (Table 7, Reaction 4, and Equation 4):
      F16BPDHAP+GAPReaction 4


      vFBA=Vf,FBA·[F16BP]Kmf,FBA+[F16BP]
      (Eq. 4)


      Triosephosphate isomerase (TPI) (Table 8, Reaction 5, and Equation 5):
      DHAPGAPReaction 5


      vTPI=Vf,TPI·[DHAP]Kmf,TPI+[DHAP]
      (Eq. 5)


      Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Table 9, Reaction 6, and Equation 6):
      GAP+NAD+Pi13BPG+NADH+H+Reaction 6


      vGAPDH=Vf,GAPDH*GAPKmf,GAPDH+GAP1+Km,NADNAD
      (Eq. 6)


      Table 4Parameters for hexokinase (HK)
      ParameterValueUnitsSource
      Rate of HKvHKmm/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.
      Vf, HK0.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.
      Kmf, HK0.0720mm
      • 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.
      Vr, HK1.06 × 10−4mm/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.
      Kmr, HK0.0420mm
      • 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.
      KHKm, 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.
      Table 5Parameters for phosphoglucose isomerase
      ParameterValueUnitsSource
      Rate of PGIvPGImm/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.
      Vf, PGI10.1mm/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.
      Kmf, PGI0.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.
      Vr, PGI9.60mm/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.
      Kmr, PGI0.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.
      Table 6Parameters for phosphofructokinase
      ParameterValueUnitsSource
      Rate of PFKvPFKmm/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.
      Vf, PFK1.33mm/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.
      Kmf, PFK0.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.
      KPFKm, 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.
      Table 7Parameters for fructose-bisphosphate aldolase
      ParameterValueUnitsSource
      Rate of FBAvFBAmm/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.
      Vf, FBA0.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.
      Kmf, FBA0.0380mm
      • 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.
      Table 8Parameters for triose-phosphate isomerase (TPI)
      ParameterValueUnitsSource
      Rate of TPIvTPImm/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.
      Vf, TPI5.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.
      Kmf, TPI1.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.
      Table 9Parameters for glyceraldehyde-3-phosphate dehydrogenase
      ParameterValueUnitsSource
      Rate of GAPDHVGAPDHmm/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.
      (1 + (Km, Pi)/([Pi]))−1 has been omitted from this reaction because the model does not fully track inorganic phosphate. Moreover, Ch'en et al. (10) illustrate how [Pi] increases, driving this term from 0.83 during normal conditions toward 1 under hypoxia.
      Vf, GAPDH5.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.
      Kmf, GAPDH0.0420mm
      • 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.
      KGAPDHm, NAD0.0580mm
      • 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.
      a (1 + (Km, Pi)/([Pi]))−1 has been omitted from this reaction because the model does not fully track inorganic phosphate. Moreover, Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ) illustrate how [Pi] increases, driving this term from 0.83 during normal conditions toward 1 under hypoxia.

      Cumulative distribution function (switch function)

      A model with reactions such as ours does not account for effects at very low concentrations of substrates (and, in reversible reactions, products). One consequence of this is the possibility of calculating negative products of a particular species, indicating that the model fails at this point. Two possible strategies can handle this: 1) stop the computation when a product approaches zero, or 2) introduce a pseudo-term to restrict this possibility. Although both slow the reaction rate to approach 0, the latter strategy allows us to observe future effects on the overall system. This approach has been used by others (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ,
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      ) via a cumulative distribution to represent the effects of a Boltzmann distribution. One may think of it in signal terms, where a reaction is gradually “switched” off as less molecules are available for reacting. Parameters for this function are given in Table 10 and Equation 7.
      s(q)=11+e[q]ab
      (Eq. 7)


      Phosphoglycerate kinase (PGK) (Table 11, Reaction 7, and Equation 8):
      13BPG+ADP+H+3PG+ATPReaction 7


      vPGK=Vf,PGK13BPG*s13BPGKmf,PGK*1+Km,ADPADP-Vr,PGK3PG*s13PGKmr,PGK*1+Km,ATPATP1+13BPGKmf,PGK+3PGKmr,PGK
      (Eq. 8)


      Phosphoglycerate mutase (Table 12, Reaction 8, and Equation 9):
      3PG2PGReaction 8


      vPGM=Vf,3PGM3PG*s13PGKmf,3PGM-Vr,3PGM2PG*s12PGKmr,3PGM1+3PGKmf,3PGM+2PGKmr,3PGM
      (Eq. 9)


      Enolase (ENO) (Table 13, Reaction 9, and Equation 10):
      2PGPEP+H2OReaction 9


      vENO=Vf,ENO2PG*s12PGKmf,ENO-Vr,ENOPEP*s1PEPKmr,ENO1+2PGKmf,ENO+PEPKmr,ENO
      (Eq. 10)


      Pyruvate kinase (PK) (Table 14, Reaction 10, and Equation 11):
      PEP+ADP+H+PYR+ATPReaction 10


      vPK=Vf,PKPEP*s1PEPKmf,PK*1+KPKm,ADPADP-Vr,PKPYR*s1PYRKmr,PK1+PEPKmf,PK+PYRKmr,PK
      (Eq. 11)


      Table 10Parameters for the cumulative distribution function
      ParameterValueUnitsSource
      Switch functions1(q)(Unitless)
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ,
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      [q]Current species valuemm
      a1 × 10−4mm
      b1 × 10−5mm
      Table 11Parameters for phosphoglycerate kinase
      ParameterValueUnitsSource
      Rate of PGKvPGKmm/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.
      Vf, PGK251mm/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.
      Kmf, PGK0.0210mm
      • 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.
      Vr,PGK16.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.
      Kmr, PGK0.510mm
      • 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.
      KPGKm, ATP0.00800mm
      • 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.
      KPGKm, 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.
      Table 12Parameters for phosphoglycerate mutase
      ParameterValueUnitsSource
      rate of 3PGMV3PGMmm/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.
      Vf, 3PGM11.2mm/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.
      Kmf, 3PGM0.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.
      Vr, 3PGM48.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.
      Kmr, 3PGM0.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.
      Table 13Parameters for enolase
      ParameterValueUnitsSource
      Rate of ENOvENOmm/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.
      Vf, ENO1.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.
      Kmf, ENO0.0450mm
      • 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.
      Vr, ENO2.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.
      Kmr, ENO0.0890mm
      • 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.
      Table 14Parameters for pyruvate kinase
      ParameterValueUnitsSource
      Rate of PKvPKmm/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.
      Vf, PK9.43mm/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.
      Kmf, PK0.110mm
      • 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.
      Vr, PK0.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.
      Kmr, PK10.0mm
      • 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.
      KPKm, ADP0.00268mm
      Kashiwaya et al. (60) give KPKm, ADP = 0.268 mm. However, other simulations (19, 61) do not account for a Km for ADP. We thus weakened the term due to [ADP] by decreasing KPKm, ADP by 2 orders of magnitude.
      a 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.
      ) give KPKm, ADP = 0.268 mm. However, other simulations (
      • Zhou L.
      • Salem J.E.
      • Saidel G.M.
      • Stanley W.C.
      • Cabrera M.E.
      Mechanistic model of cardiac energy metabolism predicts localization of glycolysis to cytosolic subdomain during ischemia.
      ,
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      ) do not account for a Km for ADP. We thus weakened the term due to [ADP] by decreasing KPKm, ADP by 2 orders of magnitude.
      The flux into the pentose phosphate pathway, away from G6P, is represented by G6P dehydrogenase. This flux has a maximum rate (0.0950 mm/s indicated 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.
      )) and is set to run at a percentage level equivalent to the percentage of maximum flux through phosphoglucose isomerase.
      Glucose-6-phosphate dehydrogenase (G6PD) (Table 15, Reaction 11, and Equation 12):
      G6PφReaction 11


      vG6PD=Vf,G6PD·vPGIvss,PGI
      (Eq. 12)


      Table 15Parameters for glucose-6-phosphate dehydrogenase
      ParameterValueUnitsSource
      Rate of G6PDvG6PDmm/s
      Data were calculated as a rate proportional to the rate of phosphoglucose isomerase.
      Vf, G6PD0.0950mm/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.
      a Data were calculated as a rate proportional to the rate of phosphoglucose isomerase.
      Glycerol-3-phosphate dehydrogenase is responsible for causing a flux of DHAP away from glycolysis. Wu et al. (
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: a phantom phenomenon.
      ) give the parameters for a random sequential bi bi equation to calculate this flux.
      Glycerol-3-phosphate dehydrogenase (G3PD) (Table 16, Reaction 12, and Equation 13):
      DHAP+NADHGly-3-P+NADReaction 12


      vG3PD=Vf,G3PD*DHAP*NADHKia,NADH*Km,DHAP+Km,DHAP*NADH+Km,NADH*DHAP+DHAP*NADH
      (Eq. 13)


      Table 16Parameters for glycerol-3-phosphate dehydrogenase
      ParameterValueUnitsSource
      Rate of G3PDvG3PDmm/s
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: a phantom phenomenon.
      Vf, G3PD20.0mm/s
      The maximum velocity was set at a value that allowed the species in glycolysis to rapidly stabilize near the values given by Kashiwaya et al. (60) when glucose and pyruvate concentrations were held constant. The steady-state ratio of G6P and glycerol 3-phosphate from Denton et al. (66) were used to calculate the steady glycerol 3-phosphate level from the G6P level reported by Kashiwaya et al. (60).
      Kia, NADH0.00107mm
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: a phantom phenomenon.
      Km, DHAP0.300mm
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: a phantom phenomenon.
      Km, NADH0.00284mm
      • Wu X.M.
      • Gutfreund H.
      • Lakatos S.
      • Chock P.B.
      Substrate channeling in glycolysis: a phantom phenomenon.
      a The maximum velocity was set at a value that allowed the species in glycolysis to rapidly stabilize near the values given 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.
      ) when glucose and pyruvate concentrations were held constant. The steady-state ratio of G6P and glycerol 3-phosphate from Denton et al. (
      • Denton R.M.
      • Yorke R.E.
      • Randle P.J.
      Measurement of concentrations of metabolites in adipose tissue and effects of insulin, alloxan-diabetes and adrenaline.
      ) were used to calculate the steady glycerol 3-phosphate level from the G6P level 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.
      ).
      Glucose-1,6-bisphosphate synthase is an enzyme that produces G1,6BP and 3PG from G1P and 1,3BPG (
      • Moxley M.A.
      • Beard D.A.
      • Bazil J.N.
      Global kinetic analysis of mammalian E3 reveals pH-dependent NAD+/NADH regulation, physiological kinetic reversibility, and catalytic optimum.
      ,
      • Wong L.J.
      • Rose I.A.
      Kinetic competence of a phosphoryl enzyme intermediate in the glucose-1,6-p2 synthase-catalyzed reaction. Purification, properties, and kinetic studies.
      ). However, to the best of our knowledge, the literature is lacking data regarding its quantitative role in specific tissues. We took steps to estimate the reaction as follows, but as this synthase becomes better understood, this model should be updated accordingly.
      Glucose-1,6-bisphosphate synthase (GBPS) (Table 17, Reaction 13, and Equation 14):
      13BPG+G1P3PG+G16BPReaction 13


      vGBPS=Vf,GBPS*13BPGKmf,GBPS*1+Km,G1PG1P-Vr,GBPS*3PGKmr,GBPS*1+Km,G16BPG16BP1+13BPGKmf,GBPS+3PGKmr,GBPS
      (Eq. 14)


      Table 17Parameters for glucose-1,6-bisphosphate synthase
      ParameterValueUnitsSource
      Rate of GBPSvGBPSmm/s
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Vf, GBPS10.0mm/s
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Kmf, GBPS0.0210mm
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Vr, GBPS6.00mm/s
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Kmr, GBPS0.510mm
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Km, G1P0.00800mm
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      Km, G16BP0.565mm
      All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (60).
      a All Km values for this reaction were treated as identical to the corresponding parameters in phosphoglycerate kinase. Because GBPS is expected to give only a minor contribution to glycolysis “parallel to” PGK, the maximum forward rate Vf, GBPS was set to approximately 1/25 of Vf, PGK. Additionally, the ratio Vf, GBPS/Vr, GBPS was approximated as Vf, 3PGM/Vr, 3PGM on account of the links between 3PGM and GBPS (
      • 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.
      ).
      In our model, G16BP is set to degrade at the same rate it is created, effectively leaving it constant.

      Glycogenesis and glycogenolysis

      Phosphoglucomutase (Table 18, Reaction 14, and Equation 15):
      G1PG6PReaction 14


      vPGM=Vf,PGM*G1PKmf,PGM-Vr,PGM*G6PKmr,PGM1+G1PKmf,PGM+G6PKmr,PGM
      (Eq. 15)


      Table 18Parameters for phosphoglucomutase
      ParameterValueUnitsSource
      Rate of PGMvPGMmm/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.
      Vf, PGM1.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.
      Kmf, PGM0.0450mm
      • 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.
      Vr, PGM1.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.
      Kmr, PGM0.670mm
      • 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.

      Inhibition of glycogenolysis and glycogenesis by AMP (Table 19)

      Although many of the parameters for glycolysis/glycogenolysis are given 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.
      ), inhibition by AMP is not included by those authors. To prevent glycogen from being employed too soon, the actual cell uses AMP to inhibit glycogen phosphorylase, which can be represented (
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      ) as shown in Equation 16.
      (1+Km,AMP[AMP])
      (Eq. 16)


      Table 19Parameters used for inhibition by AMP
      ParameterValueUnitsSource
      Km, AMP0.0160mm
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      h1.50(Unitless)
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      However, to maintain equilibrium within the glycolysis/glycogenolysis loop during perfuse conditions, our simulation required all equations to be inhibited by AMP, and so the inhibition term is included for UDP-glucose pyrophosphorylase, glycogen synthase D-form, glycogen synthase I-form, and glycogen phosphorylase. Parameters for this inhibition are given below.
      UDP-glucose pyrophosphorylase (UGP) (Table 20, Reaction 15, and Equation 17):
      G1P+UTP+H+UDPG+2·PiReaction 15


      vUGP=XUGPkf,UPG*G1P-kr,UPG*UDPG1+Km,AMPAMPh
      (Eq. 17)


      Table 20Parameters for UDP-glucose pyrophosphorylase
      ParameterValueUnitsSource
      Rate of UGPvUGPmm/s
      UGP is treated as working effectively at equilibrium in a fashion similar to creatine kinase above, where the magnitude of the relative rates were approximated by a large factor X. However, unlike creatine kinase, we did not have the equilibrium concentrations for all the species involved in UGP. Thus, we use the Keq = 4.36 given by Kashiwaya et al. (60) to approximate the relative reaction rates: Keq = [UDPG]eq[unknown]/[G1P]eq. Solving for [unknown] allowed us to obtain Keq/[unknown] = kf, UPG/kr, UPG.
      XUGP10000(Unitless)
      UGP is treated as working effectively at equilibrium in a fashion similar to creatine kinase above, where the magnitude of the relative rates were approximated by a large factor X. However, unlike creatine kinase, we did not have the equilibrium concentrations for all the species involved in UGP. Thus, we use the Keq = 4.36 given by Kashiwaya et al. (60) to approximate the relative reaction rates: Keq = [UDPG]eq[unknown]/[G1P]eq. Solving for [unknown] allowed us to obtain Keq/[unknown] = kf, UPG/kr, UPG.
      kf, UPG4.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.
      UGP is treated as working effectively at equilibrium in a fashion similar to creatine kinase above, where the magnitude of the relative rates were approximated by a large factor X. However, unlike creatine kinase, we did not have the equilibrium concentrations for all the species involved in UGP. Thus, we use the Keq = 4.36 given by Kashiwaya et al. (60) to approximate the relative reaction rates: Keq = [UDPG]eq[unknown]/[G1P]eq. Solving for [unknown] allowed us to obtain Keq/[unknown] = kf, UPG/kr, UPG.
      kr, UPG0.8801/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.
      UGP is treated as working effectively at equilibrium in a fashion similar to creatine kinase above, where the magnitude of the relative rates were approximated by a large factor X. However, unlike creatine kinase, we did not have the equilibrium concentrations for all the species involved in UGP. Thus, we use the Keq = 4.36 given by Kashiwaya et al. (60) to approximate the relative reaction rates: Keq = [UDPG]eq[unknown]/[G1P]eq. Solving for [unknown] allowed us to obtain Keq/[unknown] = kf, UPG/kr, UPG.
      Km, AMP0.0160mm
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      h1.50(Unitless)
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      a UGP is treated as working effectively at equilibrium in a fashion similar to creatine kinase above, where the magnitude of the relative rates were approximated by a large factor X. However, unlike creatine kinase, we did not have the equilibrium concentrations for all the species involved in UGP. Thus, we use the Keq = 4.36 given 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.
      ) to approximate the relative reaction rates: Keq = [UDPG]eq[unknown]/[G1P]eq. Solving for [unknown] allowed us to obtain Keq/[unknown] = kf, UPG/kr, UPG.
      Glycogen synthase D-form (GSD) (Table 21, Reaction 16, and Equation 18): 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.
      ) report the rates and constants of both the D-form and the I-form of glycogen synthase. We include both in our model.
      UDPGUDP+GlgReaction 16


      vGSD=Vf,GSD*UDPGKmf,GSD+UDPG*1+Km,AMPAMPh
      (Eq. 18)


      Table 21Parameters for glycogen synthase D-form
      ParameterValueUnitsSource
      Rate of GSDvGSDmm/s
      Data modified from Kashiwaya et al. (60).
      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.
      Kmf, 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.
      Km, AMP0.0160mm
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      h1.50(Unitless)
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      a Data modified 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.
      ).
      Glycogen synthase I-form (GSI) (Table 22, Reaction 17, and Equation 19): 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.
      ) report the rates and constants of both the D-form and the I-form of glycogen synthase. We include both in our model.
      UDPGUDP+GlgReaction 17


      vGSI=Vf,GSI*UDPGKmf,GSI+UDPG*1+Km,AMPAMPh
      (Eq. 19)


      Glycogen phosphorylase (GP) (Table 23, Reaction 18, and Equation 20):
      Glg+PiG1PReaction 18


      vGP=Vf,GP*GlgKmf,GP-Vr,GP*G1PKmr,GP1+GlgKmf,GP+G1PKmr,GP*1+Km,AMPAMPh
      (Eq. 20)


      Table 22Parameters for glycogen synthase I-form
      ParameterValueUnitsSource
      Rate of GSIvGSImm/s
      Data modified from Kashiwaya et al. (60).
      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.
      Kmf, GSI0.0800mm
      • 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.0160mm
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      h1.50(Unitless)
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      a Data modified 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.
      ).
      Table 23Parameters for glycogen phosphorylase
      ParameterValueUnitsSource
      Rate of GPvGPmm/s
      Data modified from Kashiwaya et al. (60).
      Vf, GP0.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.
      Kmf, GP0.100mm
      • 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.
      Vr, GP55.8mm/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.
      Kmr, GP5.00mm
      • 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.0160mm
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      h1.50(Unitless)
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      a Data modified 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.
      ).

      Mitochondrial oxidation of pyruvate (OxP)

      We elected to represent the mitochondrial oxidation of pyruvate as a single ensemble reaction (
      • Orlowski P.
      • Chappell M.
      • Park C.S.
      • Grau V.
      • Payne S.
      Modelling of pH dynamics in brain cells after stroke.
      ,
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      ) because we primarily emphasize glycolysis, and the reactions are more heavily represented during anaerobic considerations. Equations to directly model pyruvate dehydrogenase, the citric acid cycle, electron transport chain, and related reactions and ion exchanges are included in the estimation below. Other models (
      • 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.
      ,
      • Wu F.
      • Yang F.
      • Vinnakota K.C.
      • Beard D.A.
      Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
      ,
      • Zhou L.
      • Salem J.E.
      • Saidel G.M.
      • Stanley W.C.
      • Cabrera M.E.
      Mechanistic model of cardiac energy metabolism predicts localization of glycolysis to cytosolic subdomain during ischemia.
      ) consider the details of various aspects of these components, and future efforts may incorporate equations from those models. Parameters for this equation are given in Table 24. (See Reaction 19 and Equation 21).
      15·ADP+15·Pi+3·O2+PYR+NADH15·ATP+3·CO2+NADReaction 19


      vOxP=vPK·([ADP][ADP]+Km,ADP)·([O2][O2]+Km,O2)
      (Eq. 21)


      Table 24Parameters for the mitochondrial oxidation of pyruvate
      ParameterValueUnitsSource
      Rate of OxPvOxPmm/s
      This reaction is modified from Cloutier et al. (61). We removed their “switching function” (Equation 7) because of its handling of ADP, which is already accounted for. Additionally, because we are abbreviating reactions following pyruvate kinase, we likewise removed the term for activation by pyruvate in this equation, instead opting to match the maximum velocity to the velocity of pyruvate generation vPK.
      Km, ADP5.00mm
      • 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, O20.00297mm
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      a This reaction is modified from Cloutier et al. (
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      ). We removed their “switching function” (Equation 7) because of its handling of ADP, which is already accounted for. Additionally, because we are abbreviating reactions following pyruvate kinase, we likewise removed the term for activation by pyruvate in this equation, instead opting to match the maximum velocity to the velocity of pyruvate generation vPK.

      ATP consumption (ATPh)

      The live cardiomyocyte expresses a variable ATP consumption rate as a result of various possible feedback mechanisms. Our model explores how well the metabolic system behaves across the range of possible consumption rates by varying a set rate of hydrolysis over our simulations. We use a base ATP consumption rate as reported on beating canine heart experiments (
      • 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.
      ). Each of our simulations then explores some percentage rATP of this rate. See Reaction 20 and Equation 22.
      ATPADP+Pi+H+Reaction 20


      vATPh=rATP·VRcell,cyto·XATPh·s2(ATP)
      (Eq. 22)


      Parameters for ATP consumption are given in Table 25. The s2 cumulative distribution function merely serves to halt the fixed consumption rate once ATP falls to a physiologically low value (see the discussion of Equation 7). s2 can be considered as essentially unity, and ATP is essentially constant because our results and the analysis do not consider [ATP] less than 0.007 mm (parameters for s2 are given in Table 26). See Equation 23.
      s2(q)=11+e[q]ab
      (Eq. 23)


      Table 25Parameters for ATP consumption
      ParameterValueUnitsSource
      Rate of ATPhvATPhmm/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.
      rATP
      rATP represents the fixed percentage of ATP consumption, which is one of the two values changed from simulation to simulation.
      (Unitless)NA
      VRcell, cyto1.47Ratio of cellular to cytoplasmic volume
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      XATPh0.390mmol/(s·cellular liter)
      • 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.
      a rATP represents the fixed percentage of ATP consumption, which is one of the two values changed from simulation to simulation.
      Table 26Parameters for the switch function for ATP consumption
      ParameterValueUnitsSource
      Switch functions2(q)(Unitless)
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ,
      • Cloutier M.
      • Bolger F.B.
      • Lowry J.P.
      • Wellstead P.
      An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements.
      [q]Current species valuemm
      a1 × 10−3mm
      b1 × 10−4mm

      ATP buffering

      Two enzymes in the cytoplasm contribute to buffering ATP levels. Creatine kinase transfers a phosphate group to ADP, whereas adenylate kinase transfers a phosphate group from one ADP to another. Both reactions are considered to work at quasi-equilibrium when compared with the other reactions in the model, and so their reaction velocities are driven at some high value to achieve their equilibrium ratio (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ).
      Creatine kinase (CK) (Table 27, Reaction 21, and Equation 24):
      ADP+CrP+H+ATP+CrReaction 21


      vCK=XCK*KCK*ADP*CrP*H+-ATP*Cr
      (Eq. 24)


      Adenylate kinase (AK) (Table 28, Reaction 22, and Equation 25):
      2·ADPATP+AMPReaction 22


      vAK=XAK×(KAK·[ADP]2[AMP]·[ATP])
      (Eq. 25)


      Table 27Parameters for creatine kinase
      ParameterValueUnitsSource
      Rate of CKvCKmm/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      XCK1 × 104mm/(s·mm2)
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      KCK1.66 × 10−61/mm
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      Table 28Parameters for adenylate kinase
      ParameterValueUnitsSource
      Rate of AKvAKmm/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      XAK1 × 104mm/(s·mm2)
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      KAK1(Unitless)
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.

      Effects of ion transport on pH

      Two equations of our model, creatine kinase (Equation 24) and lactic acid association (Equation 35), depend on the number of hydrogen ions available, which requires some estimation of pH in the cell.
      The model incorporates a preliminary accounting of the effect that ion flow has on pH (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ). Cl-HCO3 exchange (AE), Cl-OH exchange (CHE), Na+-H+ exchange (NHE), and Na+-HCO3 symport (NHS) can all be represented as functions of pH; the following Equations 2629 in their long form are taken directly from Ch'en et al. (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ), with a sign change where appropriate for the direction of relevant ion flux. One caveat: the ion fluxes J have empirical formulas given to represent experimental data. The number of digits in each term do not represent the level of confidence in the accuracy of the data but are given to allow researchers to compare models avoiding numerical rounding errors for the curves.
      In addition, [HCO3] is estimated as per Lagadic-Gossmann et al. (
      • Lagadic-Gossmann D.
      • Buckler K.J.
      • Vaughan-Jones R.D.
      Role of bicarbonate in pH recovery from intracellular acidosis in the guinea-pig ventricular myocyte.
      ) as shown in Equation 30,
      JAE=160×-15.2266606471×pH4+302.2590169999×pH3-1823.1533057568×pH2+1976.4960115099×pH+8383.533719598Meq·l-1·min-1
      (Eq. 26)


      JCHE=160×-0.272561314×pH4+12.313011154×pH3-181.7704305807×pH2+1108.6191429405×pH-2422.8396631585Meq·l-1·min-1
      (Eq. 27)


      JNHE=160×20.6092567224×pH4-606.5562860276×pH3+6701.065336577×pH2-32930.5476482116×pH+60727.9345421164Meq·l-1·min-1
      (Eq. 28)


      JNHS=160×2.3290050022×pH3-45.1765173617×pH2+286.7706982101×pH-592.1682240141Meq·l-1·min-1
      (Eq. 29)


      [HCO3]=[HCO3e]×10pHpHe
      (Eq. 30)


      Using these parameters, and |[HCO3]| (the magnitude of the concentration of [HCO3] in mm), we can estimate the pH change (
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      ) as shown in Equation 31
      dpHdt=-0.8×vATPase+0.4×vCK+JAE+JCHE+JNHE+JNHS-vLHX+vLA+vLDH-vHK-vPFK-vGAPDH+vPK2.3×HCO3--28×pH+222.6
      (Eq. 31)


      where you get Equation 32,
      [H+]=1000×10pH(in mM)
      (Eq. 32)


      Lactate transport, lactate dehydrogenase, and lactic acid association

      Lactate affects both the pH and pyruvate levels. The presence of lactate evolves according to the following reactions.
      Lactate dehydrogenase (LDH) (Table 29, Reaction 23, and Equation 33):
      PYR+NADH+H+L+NADReaction 23


      vLDH=Vf,LDH*PYRKmf,LDH+PYR*1+KLDHm,NADHNADH
      (Eq. 33)


      Lactate/H+ cotransporter (LHX) (Table 30, Reaction 24, and Equation 34):
      Le+He+L+H+Reaction 24


      vLHX=Vf,LHX·[Le]Kmf,LHX+[Le]Vr,LHX·[L]Kmr,LHX+[L]
      (Eq. 34)


      Lactic acid association (LA) (Table 31, Reaction 25, and Equation 35):
      L+H+LaHReaction 25


      vLA=XLA×([L]·[H+]KLA[LaH])
      (Eq. 35)


      Table 29Parameters for lactate dehydrogenase
      ParameterValueUnitsSource
      Rate of LDHvLDHmm/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.
      Vf, LDH23.9mm/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.
      Kmf, LDH0.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.
      KLDHm, NADH0.00100mm
      • 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.
      Table 30Parameters for lactate/H+ cotransporter
      ParameterValueUnitsSource
      Rate of LHXvLHXmm/s
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Vf, LHX0.0482mm/s
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Kmf, LHX2.20mm
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Vr, LHX0.182mm/s
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Kmr, LHX6.92mm
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      Table 31Parameters for lactic acid association
      ParameterValueUnitsSource
      Rate of LAvLAmm/s
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      KLA1.259 × 10−7mm
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.
      XLA10,0001/s
      • Ch'en F.F.
      • Vaughan-Jones R.D.
      • Clarke K.
      • Noble D.
      Modelling myocardial ischaemia and reperfusion.

      Oxygen transport

      The oxygen level in the blood vessel supplying the myocyte is one of the two fixed variables in the simulation. (The second is the ATP consumption rate.) After 100 s, the blood vessel oxygen concentration is dropped from its prescribed initial value (Table 32) to any desired level of hypoxia.
      Table 32Initial oxygen value
      SpeciesInitial valueSource
      mm
      Oxygen, in blood vesselO2, v0.133
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      The partial pressure at the arterial entry is approximately P1 = 100 mm Hg (24). This is found as a concentration via O2, v, = αP1 + Hct × CHb × SHb, where SHb = P1nH/(P1nH + P50, HbnH), α = 1.3 × 10−6 m mm Hg−1, Hct = 0.45, CHb = 0.0213 mol, P50, Hb = 30 mm Hg, nH = 2.55; however, it is the dissolved oxygen that is immediately available for membrane permeation, thus [O2, v] = αP1.
      a The partial pressure at the arterial entry is approximately P1 = 100 mm Hg (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ). This is found as a concentration via O2, v, = αP1 + Hct × CHb × SHb, where SHb = P1nH/(P1nH + P50, HbnH), α = 1.3 × 10−6 m mm Hg−1, Hct = 0.45, CHb = 0.0213 mol, P50, Hb = 30 mm Hg, nH = 2.55; however, it is the dissolved oxygen that is immediately available for membrane permeation, thus [O2, v] = αP1.
      Our model accounts for membrane permeability of oxygen across the vessel-extracellular space barrier and across the cell membrane. The mitochondrion is modeled as sharing the same oxygen pool as the cytoplasm. The value of oxygen for the blood vessel compartment is held fixed and allows the oxygen to enter new compartments, each of which is approximated as being well mixed.
      Permeation of O2 (vessel-extracellular, O2ve) (Table 33, Reaction 26, and Equation 36):
      O2,vO2,eReaction 26


      vO2ve=PAveVe([O2,v][O2,e])
      (Eq. 36)


      Permeation of O2 (extracellular-cellular, O2ec) (Table 34, Reaction 27, and Equation 37):
      O2,eO2,cReaction 27


      vO2ec=PAecVc([O2,e][O2,c])
      (Eq. 37)


      Table 33Parameters for oxygen transport from blood vessel to extracellular compartment
      ParameterValueUnitsSource
      Rate of O2vevO2vemm/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      PA50.0liter/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      Ve0.241/0.06841
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      a SBML automatically handles volume sizes across compartments by calculating total amount of substance crossed, i.e. concentration times an implicit compartment volume; we have added the ratio of the compartments here for ease of reading. Beard (
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      ) gives the relative volumes of each compartment. Calculating these relative to the cytoplasm, we have Vextracellular/Vcytoplasm = 0.241 and Vvessel/Vcytoplasm = 0.0684.
      Table 34Parameters for oxygen transport from extracellular compartment to cellular compartment
      ParameterValueUnitsSource
      Rate of O2ecvO2ecmm/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      PA10.0liter/s
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      Vc
      As above, SBML automatically handles volume sizes across compartments. The relative volume Vextracellular/Vcytoplasm = 0.241 is determined from values given in Ref. 24.
      1/0.2411
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      a As above, SBML automatically handles volume sizes across compartments. The relative volume Vextracellular/Vcytoplasm = 0.241 is determined from values given in Ref.
      • Beard D.A.
      Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
      .

      Myoglobin-oxygen binding (MB) (Table 35)

      Although the dynamics for myoglobin-enhanced transport are on a fast time scale for the results in this model, they are nonetheless incorporated to facilitate future simulations (
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      ). See Reaction 28 and Equation 38.
      Mb+O2MbO2Reaction 28


      vMB=ka·[Mb]·[O2]ka·[MbO2]
      (Eq. 38)


      Table 35Parameters for myoglobin binding (MB) to oxygen
      ParameterValueUnitsSource
      Rate of MBvMBmm/s
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      ka1.54 × 1041/(mm·s)
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.
      kb60.01/s
      • Endeward V.
      The rate of the deoxygenation reaction limits myoglobin- and hemoglobin-facilitated O2 diffusion in cells.

      Author contributions

      A. D. M. wrote the paper and developed the simulations. A. D. M. and C. F. D. analyzed simulation results, reviewed the results, and approved the final version of the manuscript.

      Acknowledgments

      We thank Stephen Payne and Piotr Orlowski of the University of Oxford, Oxford, United Kingdom, and David Sosnovik of Harvard Medical School for helpful discussions. We also thank Huey Eng Chua of Nanyang Technological University in Singapore for sharing some early modeling using Cell Designer.

      References

        • United States Burden of Disease Collaborators
        JAMA. 2013; 310: 591-608
        • World Health Organization
        World Health Statistics 2014. World Health Organization, Geneva, Switzerland2014
        • Yellon D.M.
        • Hausenloy D.J.
        Myocardial reperfusion injury.
        N. Engl. J. Med. 2007; 357: 1121-1135
        • Hausenloy D.J.
        • Yellon D.M.
        Myocardial ischemia-reperfusion injury: a neglected therapeutic target.
        J. Clin. Invest. 2013; 123: 92-100
        • Kalogeris T.
        • Baines C.P.
        • Krenz M.
        • Korthuis R.J.
        Cell biology of ischemia/reperfusion injury.
        Int. Rev. Cell Mol. Biol. 2012; 298: 229-317
        • Schwartz Longacre L.
        • Kloner R.A.
        • Arai A.E.
        • Baines C.P.
        • Bolli R.
        • Braunwald E.
        • Downey J.
        • Gibbons R.J.
        • Gottlieb R.A.
        • Heusch G.
        • Jennings R.B.
        • Lefer D.J.
        • Mentzer R.M.
        • Murphy E.
        • Ovize M.
        • et al.
        New horizons in cardioprotection: recommendations from the 2010 National Heart, Lung, and Blood Institute Workshop.
        Circulation. 2011; 124: 1172-1179
        • Ferdinandy P.
        • Hausenloy D.J.
        • Heusch G.
        • Baxter G.F.
        • Schulz R.
        Interaction of risk factors, comorbidities, and comedications with ischemia/reperfusion injury and cardioprotection by preconditioning, postconditioning, and remote conditioning.
        Pharmacol. Rev. 2014; 66: 1142-1174
        • Bell R.M.
        • Bøtker H.E.
        • Carr R.D.
        • Davidson S.M.
        • Downey J.M.
        • Dutka D.P.
        • Heusch G.
        • Ibanez B.
        • Macallister R.
        • Stoppe C.
        • Ovize M.
        • Redington A.
        • Walker J.M.
        • Yellon D.M.
        9th Hatter Biannual Meeting: position document on ischaemia/reperfusion injury, conditioning and the ten commandments of cardioprotection.
        Basic Res. Cardiol. 2016; 111: 41
        • Kalogeris T.
        • Baines C.P.
        • Krenz M.
        • Korthuis R.J.
        Ischemia/reperfusion.
        Compr. Physiol. 2016; 7: 113-170
        • Ch'en F.F.
        • Vaughan-Jones R.D.
        • Clarke K.
        • Noble D.
        Modelling myocardial ischaemia and reperfusion.
        Prog. Biophys. Mol. Biol. 1998; 69: 515-538
        • Cortassa S.
        • Aon M.A.
        • O'Rourke B.
        • Jacques R.
        • Tseng H.-J.
        • Marbán E.
        • Winslow R.L.
        A computational model integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte.
        Biophys. J. 2006; 91: 1564-1589
        • 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.
        Biophys. J. 2013; 105: 1045-1056
        • Korzeniewski B.
        • Noma A.
        • Matsuoka S.
        Regulation of oxidative phosphorylation in intact mammalian heart in vivo.
        Biophys. Chem. 2005; 116: 145-157
        • Tepp K.
        • Timohhina N.
        • Chekulayev V.
        • Shevchuk I.
        • Kaambre T.
        • Saks V.
        Metabolic control analysis of integrated energy metabolism in permeabilized cardiomyocytes–experimental study.
        Acta Biochim. Pol. 2010; 57: 421-430
        • Wu F.
        • Yang F.
        • Vinnakota K.C.
        • Beard D.A.
        Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology.
        J. Biol. Chem. 2007; 282: 24525-24537
        • Yaniv Y.
        • Stanley W.C.
        • Saidel G.M.
        • Cabrera M.E.
        • Landesberg A.
        The role of Ca2+ in coupling cardiac metabolism with regulation of contraction.
        Ann. N.Y. Acad. Sci. 2008; 1123: 69-78
        • 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.
        J. Physiol. 2008; 586: 4193-4208
        • Karlstädt A.
        • Fliegner D.
        • Kararigas G.
        • Ruderisch H.S.
        • Regitz-Zagrosek V.
        • Holzhütter H.-G.
        CardioNet: a human metabolic network suited for the study of cardiomyocyte metabolism.
        BMC Syst. Biol. 2012; 6: 114
        • Zhou L.
        • Salem J.E.
        • Saidel G.M.
        • Stanley W.C.
        • Cabrera M.E.
        Mechanistic model of cardiac energy metabolism predicts localization of glycolysis to cytosolic subdomain during ischemia.
        Am. J. Physiol. Heart Circ. Physiol. 2005; 288: H2400-H2411
        • Aon M.A.
        • Cortassa S.
        Mitochondrial network energetics in the heart.
        Wiley Interdiscip. Rev. Syst. Biol. Med. 2012; 4: 599-613
        • 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.
        Ann. N.Y. Acad. Sci. 2010; 1188: 133-142
        • Lopaschuk G.D.
        • Ussher J.R.
        • Folmes C.D.
        • Jaswal J.S.
        • Stanley W.C.
        Myocardial fatty acid metabolism in health and disease.
        Physiol. Rev. 2010; 90: 207-258
        • Fillmore N.
        • Lopaschuk G.D.
        Malonyl CoA: a promising target for the treatment of cardiac disease.
        IUBMB Life. 2014; 66: 139-146
        • Beard D.A.
        Modeling of oxygen transport and cellular energetics explains observations on in vivo cardiac energy metabolism.
        PLoS Comput. Biol. 2006; 2: e107
        • Rumsey W.L.
        • Pawlowski M.
        • Lejavardi N.
        • Wilson D.F.
        Oxygen pressure distribution in the heart in vivo and evaluation of the ischemic “border zone.”.
        Am. J. Physiol. 1994; 266: H1676-H1680
        • Seiler C.
        • Stoller M.
        • Pitt B.
        • Meier P.
        The human coronary collateral circulation: development and clinical importance.
        Eur. Heart J. 2013; 34: 2674-2682
        • Das A.M.
        • Harris D.A.
        Control of mitochondrial ATP synthase in heart cells: inactive to active transitions caused by beating or positive inotropic agents.
        Cardiovasc. Res. 1990; 24: 411-417
        • Rolfe D.F.
        • Brown G.C.
        Cellular energy utilization and molecular origin of standard metabolic rate in mammals.
        Physiol. Rev. 1997; 77: 731-758
        • O'Gara P.T.
        • Kushner F.G.
        • Ascheim D.D.
        • Casey Jr., D.E.
        • Chung M.K.
        • de