Modeling oxygen requirements in ischemic cardiomyocytes

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.

Globally, ischemic heart disease was the leading cause of death in 2012 and the most rapidly increasing cause of death during 2000 -2012 (1,2). 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 (3). 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 (4 -9).
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 (10 -16). However, only Ch'en et al. (10) 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. (17) represent the onset of ischemia from normal conditions (for 30 s). An alternative approach is taken by Karlstädt et al. (18) 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. (19) 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. (19) 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, mitochon-drial 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 (20). One focal aspect of this debate is the mechanism by which ATP consumption is regulated, whether by calcium, inorganic phosphate, creatine, or otherwise (21). 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 (22). 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 (23). 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.
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 (24) 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  (22) 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.
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. (10) used a completely anaerobic model and found [ATP] to last roughly 800 s. Wu et al. (17) provide another point of comparison regarding creatine phosphate (CrP) 2 : 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. (17) are not extended to full depletion of CrP).
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. 2 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.  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 O 2 levels below a fairly sharp threshold. As one may expect, if we decrease the energy consumption rate under low O 2 conditions, [ATP] no longer plummets (Fig. 4B).
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 [O 2 ] and ATP consumption rates that lead to these [ATP] thresholds. The ATP consumption rates from these simulations and the resulting extracellular [O 2, e ] are plotted in Fig. 5.
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 O 2 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 O 2 (about 12.8 mm Hg). This falls just under the venous P O2 values of 15-20 mm Hg reported in the literature (24,25). Additionally, the fact that such small changes in oxygen can result in drastically different metabolic outcomes is consistent  Figure 5. Curves 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 Fig. 4, 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. with the well known fact that patients with significant collateral circulation upon X-ray angiography have improved outcomes following reperfusion (26). Given this map, we can begin to identify some parameters that are key for the survival of the cardiomyocyte. Das and Harris (27) 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 (27) indicate that a non-contracting cell consumes 56% of the energy that a contracting cell consumes (see Fig. 6).
Rolfe and Brown (28) 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.
Combining the ATP state information in Fig. 5 with the ATP consumption data in Figs. 6 and 7 to create Fig. 8, we find that when the extracellular O 2 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.
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 nonbeating 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

Discussion
The human heart and the cardiomyocyte are incredibly resilient. If an artery is blocked and insufficient O 2 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 O 2 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 O 2 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 (29).
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 O 2 where no change in [ATP] is detected. This means that ischemic tissue with a gradient in O 2 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 [O 2 ] gradient.
Others (30,31) 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 (31) 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 (30).
The spatial extent of cells transitioning from high to low [ATP] can occur over a range of 0.0015 mM O 2 . 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 partic-ularly susceptible to perturbations such as reperfusion. At first glance, the range of [O 2 ] 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 (32). 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 (33).
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 (26). 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 Figs. 5 and 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 (34) 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 O 2 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) (5,30,35). Recent studies (36) 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 (37). 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 (38) that there is a very substantial (45%) decrease of ATP following reperfusion.
Weakened cardiomyocytes may have less resources to handle ROS loads (39). The role of antioxidants, such as superoxide dismutase (SOD), during reperfusion remains unclear. Heart outcomes appear to improve with increased activity of SOD (40,41) as well as with administered antioxidants (42); yet SOD appears to degrade quickly when administered intravenously (43), and the level of SOD activity throughout reperfusion remains unclear (44). 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 (45).
Another important player in both ischemia and ischemia/ reperfusion injury is the calcium gradient across the inner membrane of the mitochondria. For example, Ca 2ϩ -ATPase pumps will weaken as ATP is decreased (46). 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 (47). 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 Ca 2ϩ 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 (17,18). 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 (10,12,15,48,49). One illustration of the effects due to these changing concentrations is demonstrated by Vinnakota et al. (49) 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 byproducts 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 K m, 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 (50). These include the following: elements of mitochondrial structure, such as the membrane permeability transition pore (51); ion regulation and transport, including mitochondrial potassium and calcium (52); pools of key molecular species, for example NAD/NADH (53); and information on molecular behavior, such as the details regarding ROS generation (12,54). 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 (55). For example, the mitochondrial model iAS253 for flux balance analysis by Smith and Robinson (56) was employed by Chouchani et al. (37) 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) (57,58). Reactions were represented as a system of coupled ordinary differential equations and then evaluated in the CellDesigner Simulation environment (59).

Biochemical species
The model includes 56 reactions to explicitly calculate 32 biochemical species while including the effect of other species as detailed in Tables 1-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 1 Species evaluated in the model Kashiwaya et al. (60) report the concentrations of metabolites found in working hearts from Wistar rats prepared with a modified Langendorff technique. The values from Ch'en et al. (10) listed here are from hearts prepared with a buffer solution containing 10 mM glucose. a 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. b Data were calculated from Denton et al. (66) and Kashiwaya et al. (60). c Data were calculated using the value for L Ϫ e and equations in Ch'en et al. (10). d Endeward (65) identifies ͓Mb tot ͔ ϭ 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 ͓O 2 ͔. e Early iterations of the model were run with the given ͓O 2, v ͔ value, and observed to find the steady ͓O 2, e ͔ and ͓O 2 ͔.   (15) 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 (60).

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. (60). In our model, reactions are generally represented as reversible Michaelis-Menten equations with modifiers from secondary reactants and products. Hexokinase (HK) ( Triosephosphate isomerase (TPI) ( Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) ( Table  9, Reaction 6, and Equation 6):       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 (10, 61) 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. Phosphoglycerate kinase (PGK) ( The flux into the pentose phosphate pathway, away from G6P, is represented by G6P dehydrogenase. This flux has a Glycerol-3-phosphate dehydrogenase is responsible for causing a flux of DHAP away from glycolysis. Wu et al. (62) give the parameters for a random sequential bi bi equation to calculate this flux. Glycerol-3-phosphate dehydrogenase (G3PD) ( Glucose-1,6-bisphosphate synthase is an enzyme that produces G1,6BP and 3PG from G1P and 1,3BPG (53,63). 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): (Eq. 14) In our model, G16BP is set to degrade at the same rate it is created, effectively leaving it constant.

Glycogenesis and glycogenolysis
Phosphoglucomutase (  (Table 19)-Although many of the parameters for glycolysis/ glycogenolysis are given by Kashiwaya et al. (60), inhibition by AMP is not included by those authors. To prevent glycogen from being employed too soon, the actual cell uses

Inhibition of glycogenolysis and glycogenesis by AMP
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) (

Mitochondrial oxidation of pyruvate (OxP)
We elected to represent the mitochondrial oxidation of pyruvate as a single ensemble reaction (48, 61) 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 (12,15,19) 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). 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 K eq ϭ 4.36 given by Kashiwaya et al. (60) to approximate the relative reaction rates: K eq ϭ ͓UDPG͔ eq ͓unknown͔/͓G1P͔ eq . Solving for ͓unknown͔ allowed us to obtain K eq /͓unknown͔ ϭ k f, UPG /k r, UPG .

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 (17). Each of our simulations then explores some percentage r ATP of this rate. See Reaction 20 and Equation 22.
Parameters for ATP consumption are given in Table 25. The s 2 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). s 2 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 s 2 are given in Table 26). See Equation 23.

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 (24). Creatine kinase (CK) (

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 (10). Cl Ϫ -HCO 3 Ϫ exchange (AE), Cl Ϫ -OH Ϫ exchange (CHE), Na ϩ -H ϩ exchange (NHE), and Na ϩ -HCO 3 symport (NHS) can all be represented as functions of pH; the following Equations 26 -29 in their long form are taken directly from Ch'en et al. (10), 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.          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 (24) gives the relative volumes of each compartment. Calculating these relative to the cytoplasm, we have V extracellular /V cytoplasm ϭ 0.241 and V vessel /V cytoplasm ϭ 0.0684.

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.
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 O 2 (vessel-extracellular, O 2ve ) (

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 (65). See Reaction 28 and Equation 38.