Modeling of Oxygen Transport and Cellular Energetics Explains Observations on In Vivo Cardiac Energy Metabolism

Observations on the relationship between cardiac work rate and the levels of energy metabolites adenosine triphosphate (ATP), adenosine diphosphate (ADP), and phosphocreatine (CrP) have not been satisfactorily explained by theoretical models of cardiac energy metabolism. Specifically, the in vivo stability of ATP, ADP, and CrP levels in response to changes in work and respiratory rate has eluded explanation. Here a previously developed model of mitochondrial oxidative phosphorylation, which was developed based on data obtained from isolated cardiac mitochondria, is integrated with a spatially distributed model of oxygen transport in the myocardium to analyze data obtained from several laboratories over the past two decades. The model includes the components of the respiratory chain, the F0F1-ATPase, adenine nucleotide translocase, and the mitochondrial phosphate transporter at the mitochondrial level; adenylate kinase, creatine kinase, and ATP consumption in the cytoplasm; and oxygen transport between capillaries, interstitial fluid, and cardiomyocytes. The integrated model is able to reproduce experimental observations on ATP, ADP, CrP, and inorganic phosphate levels in canine hearts over a range of workload and during coronary hypoperfusion and predicts that cytoplasmic inorganic phosphate level is a key regulator of the rate of mitochondrial respiration at workloads for which the rate of cardiac oxygen consumption is less than or equal to approximately 12 μmol per minute per gram of tissue. At work rates corresponding to oxygen consumption higher than 12 μmol min−1 g−1, model predictions deviate from the experimental data, indicating that at high work rates, additional regulatory mechanisms that are not currently incorporated into the model may be important. Nevertheless, the integrated model explains metabolite levels observed at low to moderate workloads and the changes in metabolite levels and tissue oxygenation observed during graded hypoperfusion. These findings suggest that the observed stability of energy metabolites emerges as a property of a properly constructed model of cardiac substrate transport and mitochondrial metabolism. In addition, the validated model provides quantitative predictions of changes in phosphate metabolites during cardiac ischemia.


Introduction
Over 30 years ago, Neely et al. [1] showed that adenosine triphosphate (ATP), adenosine diphosphate (ADP), and adenosine monophosphate (AMP) concentrations are maintained at essentially constant concentrations with changes in cardiac work in the isolated perfused heart. Twenty years ago, Balaban et al. [2] observed that energetic phosphate concentrations measured in vivo using NMR spectroscopy remain effectively constant over a range of cardiac workloads. These observations contradicted earlier models of the control of mitochondrial metabolism, which assumed that the rate of oxidative phosphorylation was regulated primarily by the availability of ADP [3,4]. Specifically, the cytoplasmic ADP concentration and the ratio of creatine phosphate (CrP) to ATP concentration was found to be approximately constant over a range of workload and rates of oxygen consumption in canine hearts [2,5]. To date, a credible validated biophysical model of the in vivo regulation of oxidative phosphorylation that explains the observed phenomena has not been established. Such a model would provide a theoretical basis for understanding how mitochondrial metabolism is regulated in response to changing ATP turnover rate while maintaining homeostatic concentrations of ATP, ADP, and CrP. Such a model would also form the basis of quantitative studies of the regulation of phosphate metabolites oxidative phosphorylation in the failing heart and other pathophysiological situations [6].
In this work, a detailed model of cardiac oxidative phosphorylation that was developed based on an extensive set of data obtained from isolated mitochondria [7] is integrated with an axially distributed model of oxygen transport and exchange with tissue [8]. The mitochondrial model includes the components of the respiratory chain, the F 0 F 1 -ATPase, adenine nucleotide translocase, and the mitochondrial phosphate transporter. The mitochondrial model is integrated with a model of ATP, ADP, AMP, and CrP metabolism in the cytoplasm, including the reactions of adenylate kinase, creatine kinase, and ATP consumption. Oxygen transport is governed by advection of blood in capillaries, passive permeation between blood and interstitial fluid and between interstitial fluid and cardiomyocytes, and metabolic consumption at Complex IV of the respiratory chain. An empiric relationship between cardiac perfusion and ATP consumption rate is determined by comparing model simulations to data published by Katz et al. [5]. The behavior of the resulting integrated model is compared to datasets published by Katz et al. [5], Portman et al. [9], and Zhang et al. [10], which report on in vivo changes in cardiac phosphate metabolites in response to changes in workload and to graded ischemia.

Overview of Oxygen Transport
Oxygen transport and energetic metabolism in cardiac tissue are simulated using a three-region one dimensionally distributed model [8,11], illustrated in Figure 1. The model assumes that within the capillary, interstitial space, and cellular (myocyte) space, the concentrations of oxygen and other metabolites vary primarily along the length of the capillary. Advective transport is modeled in the capillary region; the interstitial and cellular spaces are assumed to be stagnant (nonflowing). The cellular region is further subdivided into cytoplasmic and mitochondrial compartments, as described below. The full set of model variables is listed in Table 1, along with a brief description of the variables and the units associated with each variable. Note that oxygen concentrations are expressed as mass per unit volume, while the other intracellular species are expressed as mass per unit water volume in a given region. Each variable in the model is a function of distance along the capillary x and time t.
A representative model-predicted steady-state oxygen concentration profile is illustrated in Figure 1B, which plots the capillary, interstitial, and cellular PO 2 as functions of x/L, the normalized distance along the capillary. Partial pressures decrease from the arterial (x ¼ 0) to the venous (x ¼ L) end of the capillary and decrease from the capillary region to the cellular region.
To obtain the results illustrated in Figure 1B, the total pool of exchangeable phosphate in the cell (see Equation 26 in Materials and Methods) was set to 15 mM, flow was set to G ¼ 0.75 ml min À1 (g of tissue) À1 , and the rate of ATP consumption was set to J AtC ¼ 0.45 mmol s À1 (l cell) À1 . For all simulations, the input arterial PO 2 was assumed to be 100 mm Hg. For these values, the rates of oxygen consumption and oxygen extraction were found to be 5.3 lmol min À1 (g of Figure 1. Axially Distributed Blood-Tissue Exchange Model for Oxygen Transport (A) A three-region one dimensionally distributed model for oxygen transport is diagrammed. The three regions correspond to capillary, interstitium, and cell (myocyte). Blood flows in the capillary region. Oxygen is transported via passive permeation between the capillary, interstitial, and myocyte regions. (B) Predicted capillary, interstitial, and cellular PO 2 are plotted as functions of x/L, the normalized distance along the capillary. Simulation parameters are TPP ¼ 15 mM, G ¼ 0.75 ml min À1 (g of tissue) À1 , and J AtC ¼ 0.45 mmol s À1 (l cell) À1 . The input arterial PO 2 is assumed to be 100 mm Hg. The predicted rate of oxygen consumption and oxygen extraction are found to be 5.3 lmol min À1 (g of tissue) À1 and 76%, respectively. DOI: 10.1371/journal.pcbi.0020107.g001

Synopsis
To function properly over a range of work rates, the heart must maintain its metabolic energy level within a range that is narrow relative to changes in the rate of energy utilization. Decades of observations have revealed that in cardiac muscle cells, the supply of adenosine triphosphate (ATP)-the primary currency of intracellular energy transfer-is controlled to maintain intracellular concentrations of ATP and related compounds within narrow ranges. Yet the development of a mechanistic understanding of this tight control has lagged behind experimental observation. This paper introduces a computational model that links ATP synthesis in a subcellular body called the mitochondrion with ATP utilization in the cytoplasm, and reveals that the primary control mechanism operating in the system is feedback of substrate concentrations for ATP synthesis. In other words, changes in the concentrations of the products generated by the utilization of ATP in the cell (adenosine diphosphate and inorganic phosphate) effect changes in the rate at which mitochondria utilize those products to resynthesize ATP. tissue) À1 and 76%, respectively, values which are consistent with a moderate rate of cardiac work [12,13].

Relationship between Flow and ATP Consumption
To simulate the behavior of the model at different workloads, it is necessary to define a relationship between workload and coronary perfusion. In Katz et al. [5], cardiac flow and the rate of cardiac oxygen consumption, MVO 2 , are measured under different conditions. Yet in the integrated model, MVO 2 is not set as an input parameter; MVO 2 is computed as a function of the rate of ATP consumption J AtC . To determine a relationship between G and J AtC , the model was fit to data from Tables 2, 3, and 4 of Katz et al. [5], which report G and MVO 2 for paced hearts, hearts stimulated with epinephrine, and hearts stimulated with phenylephrine. For entry in these tables, the value of J AtC was varied until the model-predicted MVO 2 was equal to the given experimental measure. The resulting data are plotted in Figure 2A. The observed G and J AtC values approximately obey the linear relationship G ¼ 0.3326 þ (2.315)J AtC , where J AtC is expressed in units of mmol s À1 (l cell) À1 and G is computed in units of ml min À1 (g of tissue) À1 . The linear fit between the flow and ATP consumption rate is indicated by the straight line in the figure. This relationship is used to specify a value of G to use for simulating the system behavior at a given rate of J AtC in the following section. Figure 2B illustrates the predicted relationship between mitochondrial ATP production and oxygen consumption over the range a range of J AtC values from 0 to 1.7 mmol s À1 (l cell) À1 . The P/O ratio is the ratio of the rate of ATP synthesis by the F 0 F 1 -ATPase to the rate of consumption of oxygen atoms. At J AtC ¼ 0, the P/O ratio is 0 because no ATP is synthesized while a finite oxygen consumption flux exists to offset the rate of proton leak across the inner mitochondrial membrane. At high values of ATP consumption (and oxygen consumption), the P/O ratio approaches the theoretical maximum of 2.50 for the oxidative phosphorylation model of Beard [7].

Phosphate Metabolites as a Function of Workload
Using the developed linear relationship between workload and coronary flow, steady-state model predictions were Oxygen concentration in interstitial fluid mol (l ISF) À1 C O, 3 Oxygen concentration in myocyte mol (l cell) À1 P 1 Oxygen partial pressure in blood mm Hg P 2 Oxygen partial pressure in interstitial fluid mm Hg P 3 Oxygen partial pressure in myocyte mm Hg S Hb Hemoglobin saturation in red blood cells unitless S  calculated for a range of J AtC values from 0 to 1.7 mmol s À1 (l cell) À1 . Over this range, the predicted MVO 2 varies from 1.68 lmol min À1 (g of tissue) À1 to 15.7 lmol min À1 (g of tissue) À1 . Predictions of ADP and Pi concentrations and CrP/ATP ratio as functions of MVO 2 are plotted in Figure 3. In Figure 3, concentrations of metabolites are computed as the average value in the distributed model for a given steady-state condition.
The black curves in Figure 3A represent the relative cytoplasmic ADP level (normalized to the value at J AtC ¼ 0) as a function of MVO 2 for three different values of the total phosphate pool parameter (see Equation 26 in Materials and Methods). Also plotted are data from Figures 2A, 4A, and 6A of Katz et al. [5], which represent experimental observations for the three different stimulus protocols described above. It is apparent that within the range of observed variability, the model predictions match the experimental data for values of MVO 2 , less than approximately 12 lmol min À1 (g of tissue) À1 . At the highest simulated level of TPP (18 mM), the model predictions diverge from the observed data at MVO 2 values slightly lower than is predicted at TPP equal to 12 mM and 15 mM. However, regardless of the value of TPP, the model predicts that the relative ADP concentration is 2-fold to 3fold higher at MVO 2 ¼ 15 lmol min À1 (g of tissue) À1 than at the lower work rates. Although the amount of available data for extremely high workloads near MVO 2 ¼ 15 lmol min À1 (g of tissue) À1 is sparse, the observations indicate that even at these extreme work rates, the cellular ADP level does not significantly increase from the minimum value. Figure 3B plots the predicted CrP/ATP ratio in the cytoplasm as a function of MVO 2 , as well as the data from Figures 2C, 4C, and 6C of Katz et al. [5]. Model predictions at TPP ¼ 15 mM are most consistent with the experimental data. At this value for the total phosphate pool, the model predicts that the ratio CrP/ATP drops from 2 at low work rate to approximately 1 at the highest work rate. Over the range of MVO 2 from 2 to 6 lmol min À1 (g of tissue) À1 , the predicted CrP/ATP ranges from 1.89 to 1.76 at TPP ¼ 15 mM.
Model-predicted concentrations of cytoplasmic Pi as a function of work rate are plotted in Figure 3C. The model predicts that at TPP ¼15 mM, [Pi] c increases from 0.02 mM at 0 workload to 0.36 mM at MVO 2 ¼ 6 lmol min À1 (g of tissue) À1 to 3.1 mM at MVO 2 ¼ 15 lmol min À1 (g of tissue) À1 . Over the range for which the model best fits the observed data on ADP and CrP/ATP [MVO 2 ¼ 1.68 to 12 lmol min À1 (g of tissue) À1 ], the model predicts that Pi concentration increases from 0.02 mM to 1.6 mM. While it is difficult to absolutely quantify Pi concentration in vivo in the heart from 31 P-NMR spectroscopy, it is possible to estimate changes in the amount of Pi relative to other phosphate metabolites. Zhang et al. [14] report DP/CrP in the canine heart-the change in Pi relative to phosphocreatine-from baseline to maximally stimulated contraction to be 0.21. This value of DP/CrP corresponds to a change in Pi of approximately 2 mM, which is consistent with our model predictions.
In Figure 3D, the predicted MVO 2 is plotted as a function of absolute cytoplasmic ADP concentration. For comparison, data obtained from in vivo hearts of newborn (open squares) and adult (shaded triangles) from Portman et al. [9] are plotted in Figure 3D. Cytoplasmic ADP concentration is predicted to be approximately 0.1 to 0.15 mM in low and moderate workload conditions at the preferred value of TPP ¼ 15 mM. Figure 3E illustrates that cytoplasmic ATP concentration remains constant (with less than 10% variation) over the simulated range of work rate for all values of TPP. Thus the predicted variation in CrP/ATP ratio is due to changes in the CrP concentration, as illustrated in Figure 3F, which plots cytoplasm CrP as a function of MVO 2 .
Based on the above calculations, the parameter TPP is set to 15 mM for the remaining simulations.

Cytoplasmic versus Mitochondrial Metabolite Pools
As described above, the solid curves in Figure 3 represent the predicted levels of ADP, CrP/ATP, and Pi in the cytoplasm.
Comparing these values to the data of Katz et al. [5] is consistent with assumptions made in that paper that an  Table 2 of [5]), epinephrine protocol (Table 3 of [5]), and phenylephrine protocol (Table 4 of [5]), respectively. The solid line represents the best fit to the data, The predicted ratio between mitochondrial ATP production and rate of oxygen atom consumption (P/O ratio) is plotted for J AtC values from 0 to 1.7 mmol s À1 (l cell) À1 . At J AtC ¼ 0, the predicted MVO 2 is 1.68 lmol min À1 (g of tissue) À1 , corresponding to the oxygen consumption rate necessary to offset rate of proton leak across the inner mitochondrial membrane. DOI: 10 is based on two assumptions: 1) that the measured levels of both ATP and CrP represent concentrations in the cytoplasm and 2) that the creatine kinase (CK) reaction is maintained in equilibrium. This standard assumption is invoked, for example, in the studies of Balaban and coworkers [2,5,9]. However, Garlick and coworkers [15][16][17] find that both the cytoplasmic and mitochondrial ATP pools are visible in 31 P NMR spectroscopy.
If this were the case, then it would not be appropriate to assume that the total adenine nucleotide pool is in chemical equilibrium with the Cr and CrP, and consequently the reported relative ADP concentrations would not represent unbiased estimates of cytoplasmic ADP. The current model predicts that the total (cytoplasmic plus mitochondrial) ADP pool (plotted as dashed line in Figure 3A) is effectively constant over the full range of work rate simulated.

Predicted State Variables at Different Workloads
While the previous section illustrates the spatially averaged cardiac phosphate metabolite levels at different steady-state  [5]; data in (D) is obtained from Portman et al. [9]. In all panels, model curves are obtained by varying J AtC from 0 to 1.7 mmol s À1 (l cell) À1 and setting flow according to work rates, in this section predicted spatial profiles of model variables are explored at two different work rates. Specifically, simulations were performed at J AtC ¼ 0.45 mmol s À1 (l cell) À1 and G ¼ 0.75 ml min À1 (g of tissue) À1 (defined as moderate work rate) and at J AtC ¼ 0.90 mmol s À1 (l cell) À1 and G ¼ 1.25 ml min À1 (g of tissue) À1 (defined as high work rate). At the moderate work rate, oxygen consumption and extraction were computed to be 5.3 lmol min À1 (g of tissue) À1 and 76%; at the high work rate, oxygen consumption and extraction were computed to be 9.1 lmol min À1 (g of tissue) À1 and 78%. The predicted values of venous PO 2 are 18.5 mm Hg and 17.9 mm Hg for the moderate and high work rates, respectively.
Predicted spatial profiles of nine state variables are plotted in Figure 4, with solid lines and dashed lines representing data from the moderate and high work rates, respectively. From Figure 4A, it is apparent that at the high work rate, cellular PO 2 approaches 0 near the venous end of the capillary. At the moderate work rate, the minimum cellular PO 2 is 10.6 mm Hg; the minimum cellular PO 2 at the high work rate is 2.9 mm Hg. Although the cellular PO 2 varies along the length of the capillary, it is apparent that at a given workload, most of the state variables remain relatively constant throughout the cell region. Mitochondrial NADH and cytoplasmic ADP and Pi are predicted to be slightly elevated in the region of lowest oxygenation at the high work rate; DW is 2 mV to 3 mV lower in the region of low oxygen at the high work rate than it is for the majority of the myocyte space. The cytoplasmic ATP concentration ( Figure 4F) is approximately 98% of the total adenine nucleotide pool at moderate and high work rates. The cellular Pi concentration is predicted to increase from approximately 0.3 mM at moderate work to approximately 1 mM at high work.

Phosphate Metabolites during Graded Hypoperfusion
To explore the behavior of the integrated model during ischemia, when the tissue may become hypoxic to an extent greater than that observed in the simulations described above, it is not possible to treat the ATP consumption flux as a constant that does not depend on the available cytoplasmic ATP. To ensure that ATP concentration is not allowed to become negative, it is necessary to model the ATP consumption flux using an expression that approaches 0 as ATP concentration approaches 0. In this section, in which the high-energy phosphate metabolic response of cardiac tissue to graded hypoperfusion is simulated, the ATP consumption flux is modeled using the expression where the flux approaches the constant V AtC under normal conditions when ATP concentrations are more than 50 times higher than cytoplasm ADP concentrations. To simulate the model response to graded hypoperfusion, V AtC is set to 0.48 mmol s À1 (l cell) À1 and flow is varied from 0.025 to 1.25 ml min À1 (g of tissue) À1 , corresponding to the range of coronary blood flows imposed in the study by Zhang et al. [10]. The value V AtC ¼ 0.48 mmol s À1 (l cell) À1 was chosen to achieve reasonable agreement between the experimental data and model predictions. Over the range of myocardial flow, Zhang et al. [10] measured ATP, CrP, Pi, and oxymyoglobin saturation (S Mb ), in dog myocardium in vivo. These data are compared to steady-state model predictions in Figure 5. The raw data from Zhang et al. [10] were normalized so that the measured ATP concentration at normal flow was equal to 6.5 mM, the model-predicted ATP concentration for normoxic nonischemic conditions. The same scaling that was applied to the raw data for ATP was applied to the raw CrP and Pi data to obtain estimates of the molar concentrations of these species.
As expected, ATP, CrP, and S Mb levels drop and cytoplasm Pi increases as flow is reduced. In general, the order of magnitude and the qualitative trends in the model predictions match the experimental observations. However, the experimental observations show a smoother variation with flow in these four variables than is predicted by the model. The observed differences between model predictions and observed data of Figure 5 are expected to be in part accounted for by the fact that in the current model myocardial heterogeneities in flow and capillary length are not considered. This issue is discussed in greater detail in the Discussion section.

Control of Phosphate Metabolites
The major finding of this study is that the integrated model is able to reproduce experimental observations on ATP, ADP, CrP, and Pi over a range of workload and during coronary hypoperfusion. The integrated model explains metabolite levels observed at low to moderate workloads and the changes in metabolite levels and tissue oxygenation observed during graded hypoperfusion. The level of agreement between model predictions and experimental measurements is significant because the integrated model introduces only a single adjustable parameter. During a 10-fold increase in workload (rates of ATP and oxygen consumption), the cytoplasmic ATP, ADP, and CrP concentrations and mitochondrial NADH concentrations are predicted to change little. The cytoplasmic Pi concentration is predicted to increase from approximately 0.15 mM at an oxygen consumption of 4 lmol min À1 (g of tissue) À1 to approximately 1.0 mM at an oxygen consumption of 10 lmol min À1 (g of tissue) À1 . Over this range the model predictions closely match experimental measurements and indicate that cytoplasmic Pi level is a key regulator of the rate of mitochondrial respiration. At work rates higher than approximately 12 lmol min À1 (g of tissue) À1 , the model predictions of cytoplasmic ADP and CrP/ATP ratio are slightly elevated compared to low work states. Total cytoplasmic plus mitochondrial ADP is effectively constant over the entire observed range of work rate. In addition, the validated model provides quantitative predictions of changes in phosphate metabolites during cardiac ischemia. Comparison of model predictions to the data of Zhang et al. [10] provides independent validation of the model. Although these modeling results depend on the chosen value of the maximal rate of ATP consumption (V AtC ), the model predictions in the extreme limits of flow do not depend on the choice of V AtC . It is significant that the predicted Pi concentration of approximately 14 mM in the zero-flow limit and the predicted CrP concentration of approximately 12 mM at the high-flow limit are close to the observations of Zhang et al. [10]. Previous modeling studies of Korzeniewski et al. [18] have attributed the myocardium's ability to maintain homeostatic concentrations of ATP and CrP over a range of work rates to a ''parallel activation'' mechanism in which all respiratory complexes and the rate of substrate dehydrogenation are increased in synchrony with the rate of ATP consumption. However, no independent verification of the parallel activation hypothesis exists. In the current analysis, the biophysical model of cardiac oxidative phosphorylation and mitochondrial substrate transport was developed from an entirely independent set of data on isolated mitochondria [7]. The observed stability of energy metabolites emerges as a property of a properly constructed model of cardiac oxygen transport and mitochondrial metabolism. Thus the current analysis casts some doubt on the parallel activation hypothesis.
As concluded previously by Saks et al. [19,20], the current analysis predicts that Pi is an important feedback signal regulating respiration in cardiomyocytes. However, important differences exist between the current model and earlier models presented by the Saks group [20,21]. Specifically, the earlier models have applied an arbitrary scaling to previously developed models of mitochondrial oxidative phosphorylation, to match the observed dynamic range of oxygen consumption in isolated rat heart. The current model scales the model of isolated mitochondria based on independent estimates of cardiac mitochondrial volume density, eliminating an arbitrary fitting parameter necessary in the earlier models. Studies from the Saks laboratory have also suggested that coupling between mitochondrial creatine kinase and the adenine nucleotide translocase system on the mitochondrial inner membrane plays a role in maintaining metabolic homeostasis in the heart [20][21][22][23]. It will be valuable to incorporate the Vendelin-Saks model of the mitochondrial creatine kinase shuttle into the current model to explore its possible role in regulating respiration and metabolite concentrations.
The findings of the current study are consistent with the observations of Jeneson et al. [24,25] that, in skeletal muscle, ADP concentrations vary less than approximately 5-fold between minimum and maximal work rate. The current model predicts that ADP is 1.3 times higher at MVO 2 of 10 lmol min À1 (g of tissue) À1 than at MVO 2 of 1.7 lmol min À1 (g of tissue) À1 . At the maximal MVO 2 the predicted ADP concentration is 2.5 times greater than that at the minimum. The predicted cardiac Pi concentration is predicted to be too low to be detectable in NMR at low work rates [0.15 mM at MVO 2 of 4 lmol min À1 (g of tissue) À1 ] and to increase to the millimolar range at high work rates [1.0 mM at MVO 2 of 10 lmol min À1 (g of tissue) À1 ].

Mitochondrial P/O Ratio
The oxidative phosphorylation model predicts that the cardiac mitochondrial P/O ratio exceeds the value of 2 for MVO 2 greater than 7.8 lmol min À1 (g of tissue) À1 or approximately 200 ll O 2 min À1 (g of tissue) À1 . For lower values of oxygen consumption, the P/O ratio is less than 2, indicating that oxidative ATP synthesis is significantly less efficient at moderate work rates than at maximal work rates in the heart. This behavior is due to the fact that the predicted oxygen flux necessary to compensate for proton leak is nearly constant. Thus the relative contribution of proton leak to oxygen consumption increases with decreasing MVO 2 .
Since the oxygen flux associated with the proton leak is approximately one tenth of MVO max 2 , the maximal rate of oxygen consumption, the P/O ratio as a function of MVO 2 is approximately expressed as P/O ¼ 2.5 (MVO 2 À 0.1MVO max 2 )/ MVO 2 , where 2.5 is the theoretical maximum obtained by the model of oxidative phosphorylation in the absence of a leak current.

Future Considerations
It is likely that additional regulatory mechanisms that are not currently incorporated into the model may become important at high work rates. One putatively important signal is calcium. Cortassa et al. [26] and Jafri et al. [27][28][29] developed models of mitochondrial calcium handling and calcium-dependent activation of the TCA cycle. At high work rates, increased mitochondrial calcium may provide a feedforward signal to stimulate increased dehydrogenase activity. The current model, which uses a phenomenological function to drive the reduction of NAD þ , does not account for the possible role of calcium signaling in controlling cardiac respiration. Future work will integrate the detailed kinetics of the TCA cycle and calcium activation of the Jafri and Cortassa models into the current model to explore the impact of calcium signaling on cardiac metabolic regulation.
Previous analysis of oxygen transport in isolated perfused hearts has revealed that heterogeneities in microvascular geometry and flow are important in realistically simulating oxygen transport to cardiac tissue [30,31]. Yet the current model uses a single-capillary three-region model to simulate oxygen transport and metabolism. It is likely that the lack of heterogeneity in the model is responsible in part for the deviations between model predictions and experimental measurements of ATP, Pi, CrP, and S Mb during hypoperfusion. The predicted sharp transition from constant metabolite concentrations at higher flows to linearly decreasing concentrations at lower flows is a result of the model being composed of a single capillary. When flow is high enough to maintain normoxia, the energy metabolite concentrations are predicted to be stable. At lower flows, the myoglobin saturation is directly proportional to the fraction of the cell that is normoxic, which increases linearly with flow at a constant work rate. In the future, the current detailed model of cellular energetics will be ported to simulations in heterogeneous three-dimensional geometries, such as that described in [30].
The single-capillary model is used in the current work because the computational costs of the integrated model are substantial. Conducting such simulations based on a threedimensional model will require significant computational times, substantial optimization of the simulation codes, and/or increases in available computing resources.
The current model analysis of the data on flow and oxygen consumption reveals an approximately linear relationship between whole-heart cardiac perfusion and the rate of ATP consumption (illustrated in Figure 2A). Yet it has been observed that local myocardial blood flow can vary regionally by up to 10-fold in vivo [32][33][34][35]. Since oxygen [35] and fatty acid [36] uptake are tightly proportional to regional flow in the myocardium, it is proposed that ATP hydrolysis varies regionally as well. Future studies that link coronary blood flow, oxygen and substrate transport, cellular and wholeheart mechanics, and energy metabolism are expected to shed light on the physiological basis for regional heterogeneities in substrate delivery in the heart.

Summary
In sum, the integrated model of cardiac oxygen transport and mitochondrial oxidative phosphorylation is able to explain a set of data that has persisted unexplained in the literature for two decades. The model analysis predicts that Pi is a key regulator of oxidative phosphorylation in the heart, allowing ATP and CrP levels to remain relatively constant over a large range of cardiac work rate. The integrated model provides the basic framework for simulating the metabolic response of the myocardium to ischemia and hypoxia.

Materials and Methods
Governing equations for oxygen transport. Oxygen transport in the capillary is governed by advection of the blood, oxygenhemoglogin binding and unbinding, and passive permeation between the blood and the interstitial space. The governing equation for oxygen in the capillary is where C O,1 (x,t) is the total oxygen concentration in the blood (free plus hemoglobin-bound oxygen), G is blood flow to the tissue in volume per unit time per mass of tissue, q is the tissue density in mass per unit volume, L is the length of the capillary, V 1 is the volume of the capillary region, a 1 is the oxygen solubility coefficient in blood, PS 12 is the permeability of the surface area of capillary wall, and P 1 and P 2 are the oxygen partial pressures in blood and interstitial fluid, respectively. The capillary oxygen concentration C O,1 is a function of the distance along the capillary, x, and time, t. The total oxygen concentration of the capillary region is related to the partial pressure by where Hct is the hematocrit, C Hb is the concentration of oxygen binding sites in red blood cells, and S Hb is the hemoglobin saturation in red blood cells. The hemoglobin saturation curve is assumed to be governed by a Hill equation, where P 50,Hb is half the partial pressure of O 2 saturation in hemoglobin and n H is the Hill exponent. Oxygen transport in the interstitial region is governed by passive permeation of oxygen between the blood and the interstitial fluid and between the interstitial fluid and the myocyte where C O,2 (x,t) is the oxygen concentration in interstitial region, V 2 is the volume of the interstitial region, PS 23 is the permeability surfacearea product for passive permeation between interstitium region and the parenchymal cell region, and P 3 is the partial pressure of the oxygen in the cell. Oxygen concentration and partial pressure in the interstitial fluid are related by C O,2 ¼ a 2 P2, where a 2 is the oxygen solubility coefficient in the interstitial region.
In the cellular region, oxygen transport is governed by passive permeation and consumption of oxygen in mitochondria where C O,3 is the oxygen concentration in the cellular region, V 3 is the volume of the cellular region, R MitoCell is the ratio of mitochondrial volume to total cell volume, and J C4 is the flux through Complex IV in the respiratory chain in units of mass per time per unit mitochondrial volume. The factor of ½ in the second term on the RHS of Equation 6 accounts for the fact that J C4 in the mitochondrial model (described below) represents the flux of electron pairs through the respiratory chain: one O 2 molecule is consumed for each two pairs of electrons transferred through the system. Equation 6 assumes that mitochondria are uniformly and homogeneously distributed in the cellular space.
Total oxygen in the cell is the sum of free oxygen plus myoglobinbound oxygen Where C Mb is the myoglobin concentration in the cell and a 3 is the oxygen solubility in the myocyte. The equilibrium myoglobin saturation S Mb is expressed where P 50,Mb is the half-saturation partial pressure for oxygenmyoglobin binding.
Governing equations for cellular energetics. The differential equations for species other than oxygen are PLoS Computational Biology | www.ploscompbiol.org September 2006 | Volume 2 | Issue 9 | e107 In the above set of equations, the subscripts x, i, and c denote mitochondrial matrix, intermembrane (IM) space, and cytoplasm, respectively. All of the variables in this set of equations are defined in Table 1.
In addition to the state variables treated in Equations 1, 4, 5, and 8, the concentrations of several species are computed where NAD tot , Q tot , cytC tot , and A tot are the total concentrations of NAD(H), ubiquinol, cytochrome c, and adenine nucleotide in the matrix, respectively, and CR tot is the total creatine plus creatine phosphate concentration in the cytoplasm. Parameters that appear in the above equations are described in detail below. The fluxes that appear on the right-hand side of the governing equations are tabulated in Table 2. For mitochondrial species, the governing equations follow from [7]. For cytoplasmic species, the reactions modeled are ATP consumption, creatine kinase reaction, adenylate kinase reaction, and transport between the cytoplasm and the mitochondrial IM space.
Mathematical expressions for mitochondrial fluxes. Flux expressions for the mitochondrial model (based on a previously developed model [7]) are listed below. Definitions of the variables and parameters that appear in the following expressions are listed in Tables 2 and 3. Dehydrogenase flux: Complex I flux: where Complex III flux: Complex IV flux:    [37]. Values of the mitochondrial model parameters are estimated and reported in [7]. Oxidative phosphorylation model parameters have been adjusted to account for modifications from [7], as indicated in the table. Since the F 0 F 1 -ATPase reaction is maintained near chemical equilibrium in the model parameterization that was published in [7], here the F 0 F 1 -ATPase activity X F1 is set to the arbitrarily high value of 1,000 mol s À1 M À1 (l mito) À1 . In addition, the updates to the functional forms for the Complex I and III fluxes require estimation of values for the activities of X C1 and X C3 . To account for these changes, the full set of mitochondrial model parameters has been refined by repeating the fits to data from isolated mitochondrial function of Beard [7]. The agreement between the data of Bose et al. [41] and the updated model is equivalent to that reported for the previous version of the model [7].
All values for concentrations of pooled metabolites are set according to values reported in previous studies, with the exception of the total pool of exchangeable phosphate (TPP), which is estimated in the Results section. Binding constants are obtained from the literature; enzyme activities for reactions maintained near equilibrium are set to arbitrarily high values.
Computational methods. The model was coded using MATLAB; multiple steady-state solutions were obtained using the Matlab Distributed Computing Toolbox on 25-node cluster (50-CPU AMD Opteron-based server HP DL145 G2 with dual 2.6-GHz processors). Obtaining a single steady-state solution requires approximately 30 min of computing time on a single processor. A total of 150 steadystate solutions are computed to construct Figure 4, requiring less than 2 h to complete using the cluster.