Computer Modeling of Mitochondrial Tricarboxylic Acid Cycle, Oxidative Phosphorylation, Metabolite Transport, and Electrophysiology

,

A computational model of mitochondrial metabolism and electrophysiology is introduced and applied to analysis of data from isolated cardiac mitochondria and data on phosphate metabolites in striated muscle in vivo. This model is constructed based on detailed kinetics and thermodynamically balanced reaction mechanisms and a strict accounting of rapidly equilibrating biochemical species. Since building such a model requires introducing a large number of adjustable kinetic parameters, a correspondingly large amount of independent data from isolated mitochondria respiring on different substrates and subject to a variety of protocols is used to parameterize the model and ensure that it is challenged by a wide range of data corresponding to diverse conditions. The developed model is further validated by both in vitro data on isolated cardiac mitochondria and in vivo experimental measurements on human skeletal muscle. The validated model is used to predict the roles of NAD and ADP in regulating the tricarboxylic acid cycle dehydrogenase fluxes, demonstrating that NAD is the more important regulator. Further model predictions reveal that a decrease of cytosolic pH value results in decreases in mitochondrial membrane potential and a corresponding drop in the ability of the mitochondria to synthesize ATP at the hydrolysis potential required for cellular function.
Mitochondrial energy metabolism centers on the tricarboxylic acid cycle reactions, oxidative phosphorylation, and associated transport reactions. It is a system in which biochemical reactions are coupled to membrane electrophysiology, nearly every intermediate acts as an allosteric regulator of several enzymes in the system, and nearly all intermediates are transported into and out of the mitochondria via a host of electroneutral and electrogenic exchangers and cotransporters. Thus, it is a system with a level of complexity that begs for computational modeling to aid in analysis of experimental data and development and testing of quantitative hypotheses. In addition, as is demonstrated in this work, computer modeling of mitochondrial function facilitates the translation of observations made in one experimental regime (isolated ex vivo mitochondria) to another (in vivo cellular energy metabolism).
In addition, the developed model provides the basis for examining how mitochondrial energetics is controlled in vivo. The model is used to predict the roles of NAD and ADP on regulating tricarboxylic acid cycle dehydrogenase fluxes, demonstrating that NAD is the more important regulator. The mitochondrial redox state in turn is affected by cytoplasmic P i concentration, since inorganic phosphate is a co-factor for transport of tricarboxylic acid cycle substrates a substrate for the tricarboxylic acid cycle. Since ADP and P i are the biochemical substrates for oxidative phosphorylation, the primary mechanism of control of mitochondrial energy metabolism (tricarboxylic acid cycle and oxidative phosphorylation) in striated muscle is feedback of the products of ATP hydrolysis.
In a series of works, Kohn et al. Intact functional mitochondrial isolated from tissue provides a powerful tool for developing a mechanistic model of mitochondrial function. Suspensions of isolated mitochondria are readily subjected to state perturbations by introducing different substrates; mitochondrial state variables (membrane potential, redox state of cofactors, respiration rate, and intermediate concentrations) are available to measurement by a variety of methods. Therefore, a tremendous amount of valuable data shedding light on mitochondrial function is available. The challenge in making use of the available data is in building and validating a single model that is comprehensive enough to simulate a number of different experimental protocols and measurements.
This challenge is met here by building a detailed model for all of the processes illustrated in Fig. 1. This model accounts for detailed biochemical thermodynamics using a simulation approach for biochemical kinetics adapted from Vinnakota et al. (4). The electrophysiology of the mitochondrial inner membrane is based on previous models of Beard (5) and Wu et al. (6), accounting for ionic species that are transported across the membrane (such as H ϩ , ATP 4Ϫ , H 2 PO 4 Ϫ ) as variables. Although carefully estimated values for a great number of model parameters, including basic thermodynamic data of metal ion-biochemical species dissociation constants, are available from existing data bases, a large number of kinetic constants in the model have to be estimated from comparing model simulations with experimental data. In fact, the current work is unprecedented in terms of the complexity of the model and large amount of data used to build it. A total of 31 parameters values are estimated based on comparison with 25 data curves measured in isolated cardiac mitochondria from two different laboratories. Following parameterization, the developed model is validated based on comparison with additional data from both ex vivo isolated mitochondria and in vivo 31 P NMR data from skeletal muscle. The validated model is then used to investigate the major factors controlling in vivo tricarboxylic acid cycle flux and striated muscle energy metabolism.

MATERIALS AND METHODS
Here the basic approach to computer simulation and parameterization of the system illustrated in Fig. 1  dria, includes reactions occurring in three compartments: mitochondrial matrix, mitochondrial intermembrane space, and buffer space. The model incorporates tricarboxylic acid cycle fluxes, mitochondrial oxidative phosphorylation fluxes, substrate and cation transport fluxes, passive permeation fluxes, and buffer reaction fluxes. Reaction fluxes are modeled based on detailed kinetic mechanisms provided in the work of Kohn et al. (1)(2)(3)(7)(8)(9)(10)(11) and other sources and the oxidative phosphorylation model published previously (5). In total, 42 flux expressions are included in the model, including 11 tricarboxylic acid cycle fluxes, four oxidative phosphorylation fluxes, 12 substrate and cation transport fluxes across the inner mitochondrial membrane, one mitochondrial intermembrane space reaction, 11 substrate passive permeation fluxes across the outer mitochondrial membrane, and four external space reaction fluxes. All reference reactions of these fluxes, except the passive permeation fluxes, are listed in Table 1. When a reaction involves species in multiple compartments, the identifiers x, i, and c in parentheses are used to denote matrix, intermembrane, and external (buffer or cytoplasm) compartments for the species.
The model is used to simulate the kinetics of 63 state variables for the isolated mitochondrial experiments described below. These state variables include 59 biochemical reactant concentrations, matrix-free Mg 2ϩ and K ϩ , matrix pH, and inner membrane electrostatic potential. The 59 biochemical concentration variables represent 34 biochemical reactants, some of which are distributed in multiple compartments. Each of the 34 biochemical reactants is made up of several rapidly interconverting species. For example the reactant ATP is considered to be made up of the four ionic species, ATP 4Ϫ , HATP 3Ϫ , MgATP 2Ϫ , and KATP 3Ϫ . For all reactants, binding of H ϩ , K ϩ , and Mg 2ϩ is explicitly accounted for where such binding is significant.
Tricarboxylic Acid Cycle and Related Reactions-The matrix biochemical reactions considered include the tricarboxylic acid cycle reactions plus pyruvate dehydrogenase (Fig. 1, reaction 1), nucleoside diphosphokinase (reaction 10), and glutamate oxaloacetate transaminase (reaction 11). The reference chemical reactions for these biochemical reactions are tabulated in Table  1. Inhibitors and activators for the enzyme catalyzing these reactions that are considered in the model are listed in Table 2. The mechanism and corresponding mathematical expressions for each enzyme are provided in the supplemental materials, along with the kinetic parameter values used.
With the exception of enzyme activity values, all parameter values are obtained from the literature, including the modeling papers of Kohn et al. (1)(2)(3) and original experimentally based publications. In each case, our particular choice of parameter values is justified, as described in the supplemental materials. However, for many parameters, a wide range of possible values are available. We use the supplementary materials to tabulate and document alternative parameter values of kinetic constants for the tricarboxylic acid cycle enzymes. Enzyme activity values are treated as adjustable parameters, which are estimated as described below.
Enzyme regulation by metabolic intermediates is incorporated into the modeled enzyme mechanisms. Major sites of regulation in the tricarboxylic acid cycle are citrate synthase, isocitrate dehydrogenase, and ␣-ketoglutarate dehydrogenase reactions that tend to be maintained far from equilibrium (12,13). These three enzymes are strongly inhibited by product accumulation. For example, citrate acts as a competitive inhibitor against oxaloacetate, and CoA-SH acts as an uncompetitive

List of biochemical reactions
When a reaction involves species in multiple compartments, the identifiers x, i, and c are used to denote matrix, intermembrane, and external (buffer or cytoplasm) compartments for the species.

Reaction number
Enzyme Reference reaction Pyruvate-H ϩ co-transporter Hexokinase inhibitor against acetyl-CoA for citrate synthase. Both the isocitrate dehydrogenase and ␣-ketoglutarate dehydrogenase are regulated by ATP and ADP, which act as an inhibitor and an activator, respectively, at regulatory sites on these enzymes. Isocitrate dehydrogenase is an allosteric enzyme and inhibited by NADH competing against NAD; ␣-ketoglutarate dehydrogenase is inhibited by succinyl-CoA and NADH.
Oxidative Phosphorylation-The oxidative phosphorylation model of Beard and co-workers (6,14,15), including complex I, complex III, complex IV, and F o F 1 -ATPase, is extended to include succinate dehydrogenase. In previous applications using this model, the tricarboxylic acid cycle is not explicitly simulated. Instead, a phenomenological driving function is used to generate NADH and drive the electron transport system. Here, in replacing the phenomenological driving function with the biochemical model of the tricarboxylic acid cycle, the oxidative phosphorylation model parameter values are adjusted to ensure that the integrated model simultaneously matches the original data set used to identify the oxidative phosphorylation component and the kinetic data used to identify the tricarboxylic acid cycle model.
Substrate and Cation Transport Fluxes-Since the mitochondrial inner membrane is permeable to few metabolites or ions, almost all substrate and ion transport is catalyzed by specific transporters. An exception is the proton leak, which is driven by electrical potential across the inner mitochondrial membrane and is modeled here as coupled diffusion and drift using the Goldman-Hodgkin-Katz equation (16,17). In contrast to the inner membrane, the outer member is permeable to almost all small molecules and ions (13,18). All fluxes across the outer membrane are modeled as passive permeation driven by concentration gradients.
External Space Fluxes-Analysis of available data requires simulation of a variety of conditions and experimental protocols. These protocols require simulation of biochemical reactions in the external space. For the experiments of LaNoue et al. (23) (described under "Results"), hexokinase and glucose are added in the buffer medium to consume ATP and maintain mitochondrial ATP synthesis. Thus for these experiments, the hexokinase reaction is simulated. For simulations of in vivo mitochondrial energetics, creatine kinase, adenylate kinase, and ATP hydrolysis in the cytoplasm are simulated as described under "Validation." Simulation Method-The system is simulated using an approach that formally treats biochemical reactants as sums of distinct species formed by different hydrogen and metal ion binding states (4,24). The method accounts for pH and ionic dependence on enzyme kinetics and apparent equilibrium and thermodynamic driving forces for biochemical reactions. Analytical expressions for enzyme and transporter fluxes are derived based on their detailed kinetic mechanisms. Enzyme kinetics parameters, such as Michaelis-Menten constants, inhibition constants, and activation constants, are mined from a variety of sources. The detailed model equations are developed in the supplemental materials.
Parameterization Approach-All enzyme and transporter activities are treated as adjustable parameters, with values estimated based on comparison with experimental data. In total, 35 parameter values are estimated. This large number of parameters requires a large amount of relevant data for effective identification. Here, we make use of independent data sets published by LaNoue et al. (23) and Bose et al. (25). The LaNoue et al. (23) data were measured from isolated rat heart mitochondria in both resting state (state 2) and the active state (state 3), with pyruvate and malate or only pyruvate as substrates. The Bose et al. (25) data were measured from isolated pig heart mitochondria in state 2 and state 3, with glutamate and malate as substrates. The different substrates and protocols used in these experiments ensure that 2 The abbreviation used is: ANT, adenine nucleotide translocase.

Enzymes Inhibitors a Activators
Pyruvate dehydrogenase a Letters in parentheses denote types of inhibitors (C, competitive; NC, noncompetitive; UC, uncompetitive). the model is challenged by a wide range of data corresponding to diverse conditions. The total experimental data used for model identification provide 25 data curves, including 15 time courses (measure of a variable as a function of time) and 10 steady-state data sets (measure of one steady-state variable as a function of another). Based on these 25 data curves, a reasonable identification of the model is possible. Parameter values are estimated based on a Monte Carlo algorithm used to minimize the difference between model simulations and experimental data.

Experiments 1 and 2: LaNoue et al. (23) State-2 and State-3 Time Courses with Pyruvate and Malate as Substrates-La-
Noue et al. (23) used rapid quenching to measure the time courses of tricarboxylic acid cycle intermediates following the addition of malate and/or pyruvate to suspensions of isolated rat cardiac mitochondria (23). Data were reported for state 2 (with no ADP available for ATP synthesis; experiments 1 and 3) and for state 3 (with active ATP synthesis; experiments 2 and 4). State 3 was maintained in experiments 2 and 4 by including glucose and hexokinase in the buffer medium. To simulate these experiments, our mitochondrial model is initially incubated with no carbon substrate, to generate a fully oxidized initial state. At time t ϭ 0, pyruvate and malate or only pyruvate are added to the buffer, and the system is simulated for 500 or 250 s for comparison with the data of LaNoue et al. (23). Fig. 2, A-E, illustrates the state-2 time courses of pyruvate, citrate, ␣-ketoglutarate, succinate, fumarate, and malate, following the addition of 2 mM pyruvate and 5 mM malate to a suspension containing 3.4 mg of mitochondrial protein/ml of buffer. In addition, the buffer initially contains 20 mM inorganic phosphate, 5 mM magnesium ion, and 150 mM potassium ion at pH 7.2. Model simulations of the experiments are illustrated as solid lines in Fig. 2. For most reactants, the total concentration (in buffer plus mitochondria) in the suspension is provided. For citrate, time courses of both total and extramitochondrial concentrations are provided. In this experiment, pyruvate is consumed at a steady rate, and its overall concentration decreases linearly. As substrates are consumed, other intermediates build up in the system.
The state-3 experiment of LaNoue et al. (23) is identical to the state-2 experiment with the exception that 0.5 mM ADP and 40 units of hexokinase are added along with substrates pyruvate and malate at time t ϭ 0. In addition, based on the total consumed pyruvate reported by LaNoue et al. (23), we use an initial concentration of 2.5 mM in the computer simulations. Experimental measures and model simulations of the resulting time course data are illustrated in Fig. 3, A-E. In this case, pyruvate is consumed much more rapidly than in the state-2 case. The simulation predicts that the pyruvate is consumed almost linearly during the first 2 min, but the consumption rate is slowed at the end of the experiment due to limitation of available inorganic phosphate in the matrix. In addition, malate is significantly consumed in the simulation, in contrast  (23) conducted other experiments to determine the kinetics of aspartate and glutamate in resting and active states. Details for these experiments are similar to those for experiments 1 and 2 described above, with the exception that reactions are initiated by adding 1 mM pyruvate into the buffer in the absence of malate, and the system is simulated for 250 s. To simulate the experiments, aspartate is loaded into mitochondria to match the initial conditions reported by LaNoue et al. (23). Model predictions are plotted along with the experimental data points in Figs. 2F and 3F; total glutamate and aspartate concentrations in state 2 are presented in Fig. 2F, and concentrations in state 3 are plotted in Fig. 3F. The simulations show that the mitochondria convert aspartate into glutamate via glutamate oxaloace-tate transaminase following the addition of pyruvate. Since the ASP Ϫ /HGLU 0 exchanger is electrogenic, aspartate is driven out of the matrix in energized mitochondria. Simulations predict that glutamate and malate approach steady state quickly in both state-2 and state-3 simulations. The transient in state 2 is shorter than that in state 3 because of higher ASP Ϫ /HGLU 0 flux driven by higher membrane potential in state 2. The match between the data and the simulation is better for the state-3 experiment than for the state-2 experiment. Note that in the model, total glutamate plus aspartate is a conserved pool, where the total measured glutamate plus aspartate does not remain exactly constant. This inconsistency between the simulation and experiment may be a consequence of additional reactions involving glutamate and/or aspartate that are not included in the model. Specifically, the early depletion of the glutamate-aspartate pool may be due in part to deamination of glutamate by glutamate dehydrogenase, a reaction not included in the model.   Fig. 4, with modelsimulated variables plotted as solid lines. Since phosphate is required to transport malate into the matrix and is a substrate for succinyl-CoA synthetase, mitochondrial NADH and ⌬⌿ increase with [P i ] c . In fact, at [P i ] c ϭ 0, the computer model predicts that [NADH] x and ⌬⌿ ϭ 0, because no tricarboxylic acid cycle reaction flux is possible without phosphate. It is likely that some finite contamination of phosphate is present in the experiments of Bose et al. (25), even in base-line conditions with no phosphate added (25).
Simulations of the state-3 data follow the same protocol as above with the addition that ADP is added to the buffer at a concentration of 1.3 mM at t ϭ 60 s. Thus, in this case, mitochondria actively synthesize ATP when inorganic phosphate is present. As a result, respiration rates are higher, whereas mem-brane potential and NADH are lower than in state 2. The state-3 experimental data are plotted as triangles in Fig. 4, with model simulations plotted as dashed lines.
The data set illustrated in Fig. 4 was used to develop and parameterize our oxidative phosphorylation model, using a phenomenological driving function to generate NADH. When the phenomenological function is replaced by a detailed model of the tricarboxylic acid cycle, the model predictions are qualitatively similar to the original formulation of the model (5), with slightly improved predictions for the membrane potential data.
Parameter Estimation and Sensitivity-Values for the 31 adjustable parameters, listed in Table 3 Table 3 is used to perform all of the simulations illustrated in these figures.
To estimate the sensitivity to finite changes in parameter values, the sensitivity to each parameter was computed as the relative change in mean squared error due to a 10% change in a given parameter value as follows, where E* represents the minimum mean squared difference between model simulations and experimental data, and x i is the optimal value of the ith parameter. The term E i *(x i Ϯ 0.1x i ) is the error computed from setting parameter x i to 10% above and below its optimal value. The relative sensitivities to the adjustable parameters are listed in Table 3. These sensitivity values represent a measure of the degree to which the curves plotted in Figs. 2-4 are sensitive to the value of the individual parameters. A high sensitivity value indicates that changing a given parameter results in significant changes to the simulated curves used to identify the set of adjustable parameter values.
Of the 31 adjustable parameters, six are found to have sensitivities to the data of less than 1%. As indicated in the table, these low sensitivity parameters correspond to activities of enzymes and transporters that are determined to operate near equilibrium. These activities are estimated to be high enough to maintain the reactions near equilibrium, and small changes in parameter value are not expected to have a significant impact on model predictions.

Validation
To validate the model behavior, we compare its predictions with data that were not used to parameterize it. Here, we examine additional data from LaNoue et al. (23) on ␣-ketoglutarate as a function of malate concentration in the buffer and aspartate transport in state 2 and 3 (26). We also examine model predictions of in vivo ADP and P i concentrations in human skeletal muscle. Fig. 3C, ␣-ketoglutarate concentration tends to level off to an approximate constant level between 300 and 500 s in the state-3 experiment of LaNoue et al. (23). LaNoue et al. (23) reported the ␣-ketoglutarate concentration obtained with different concentrations of malate in the buffer. In Fig. 5, we compare the measured data (circles) with the model predictions (solid curve). With the exception that the initial malate concentration is varied, the simulation protocol is the same as that of Experiment 2. Model predictions are shown for the maximal ␣-ketoglutarate obtained during the simulation (t ϭ 500 s). The accumulation of ␣-ketoglutarate in the buffer is proportional to initial buffer malate concentrations and matches the experimental data points fairly well.

Relationship between Maximal [AKG] c and Initial [MAL]c in State 3-As is apparent from
In these experiments, malate in buffer enters the matrix via AKG 2Ϫ /MAL 2Ϫ , SUC 2Ϫ /MAL 2Ϫ , and MAL 2Ϫ /P i 2Ϫ exchang-ers, resulting in an increase of the matrix malate concentration. Increases in malate in the matrix accelerate tricarboxylic acid cycle turnover and ␣-ketoglutarate production rate. The simulated time courses (details not shown) show that during the initial period after malate is added into buffer, isocitrate dehydrogenase flux increases more quickly than ␣-ketoglutarate dehydrogenase flux does, resulting in an initial increase in matrix ␣-ketoglutarate. Due to the high activity of the AKG 2Ϫ /  MAL 2Ϫ exchanger, most of the produced ␣-ketoglutarate is transported into the buffer space, contributing to the overall ␣-ketoglutarate accumulation.
Aspartate Transport under Different Mitochondrial Energy States-The electrogenic ASP Ϫ /HGLU 0 exchanger is driven by both substrate concentration gradients and membrane potential (21,27,28). As shown in Figs. 2F and 3F, different steady states and time profiles of aspartate and glutamate are predicted for the state-2 and state-3 experiments. In a separate experiment, LaNoue et al. (26) investigated aspartate transport in different energy states in rat heart mitochondria. To simulate this experiment, we incubate the mitochondrial model to obtain a fully oxidized initial state as described above for Experiment 1. Buffer conditions are the same as those of Experiment 1 except that instead of pyruvate and malate, 20 mM glutamate and 1 mM malate are added into the buffer at time t ϭ 0. After 60 s, 0.1 mM ADP and 5 units of hexokinase are added into the buffer to activate state 3. The simulations are terminated at time t ϭ 150 s. Fig. 6A Fig. 6B shows the model-predicted membrane potential during this experiment. The predicted final matrix aspartate concentration is 0.24 nmol/mg protein, whereas the corresponding measured value is ϳ0.37 nmol/mg protein. However, over the course of the experiment, the experimentally measured matrix aspartate is significantly higher than the model prediction, perhaps in part due to binding of aspartate to matrix proteins not accounted for in the model. The simulations predict that the aspartate efflux rate is higher in state 3 than in state 2. When ADP and hexokinase are added into the buffer, state 3 is activated with membrane potential lowered from ϳ200 to ϳ180 mV. The decrease in membrane potential results in decreased aspartate transport out of the inner membrane. Meanwhile, the glutamate oxaloacetate transaminase reaction (proceeding in the reverse direction) is accelerated due to decreased levels of reducing equivalents in the matrix, resulting in higher production rate of the matrix aspartate. Consequently, the matrix aspartate concentration is elevated, leading to higher aspartate efflux rate compared with state 2, even with a somewhat lower membrane potential.
In Vivo Concentrations of Phosphate Metabolites in Skeletal Muscle-We have demonstrated that our model of oxidative phosphorylation, integrated into a model of cellular energetics, mimics the observed relationship between work rate (rate of oxygen consumption or rate of ATP hydrolysis) and ADP, P i , and phosphocreatine measured using 31 P NMR in cardiac (29) and skeletal muscle (6,29) in vivo. Here we show that the current model remains capable of explaining the observed data as well as or better than the simpler model. To simulate oxidative metabolism in vivo, the cytoplasmic ATP hydrolysis, creatine kinase, and adenylate kinase reactions are included in the model as described previously (6). To supply the tricarboxylic acid cycle, the cytoplasmic pyruvate concentration is held fixed at 0.06 mM.
Model predictions of the relationship between phosphate metabolites and ATP hydrolysis rate are plotted as curve 1 in Fig. 7, along with in vivo 31 P NMR spectroscopy data collected from exercising human flexor forearm muscle in healthy subjects (30). The current model prediction matches the experimental data better than simulation results of our previous skeletal muscle model (6), particularly at the lowest work rates. These findings, consistent with our previous studies, demonstrate that the observed relationships between workload and phosphate metabolites in skeletal muscle are explained by a model in which ATP synthesis is primarily controlled by feedback of substrate (ADP and inorganic phosphate) concentrations.

Predictions
Regulation of Tricarboxylic Acid Cycle Fluxes by NAD and ADP-Although a large number of regulatory mechanisms are simulated in the model, the primary control of tricarboxylic acid cycle fluxes is expected to be through cellular phosphorylation potential and redox state. There can be no net flux through the tricarboxylic acid cycle when concentration of either NAD or ADP, which serve as substrates for reactions in the cycle, is zero. Thus, when the ratios [ATP]/[ADP] and [NADH]/[NAD] are high, we expect the tricarboxylic acid cycle reaction fluxes to be inhibited by simple mass action. In addition, the allosteric inhibition of several enzymes (e.g. inhibition of pyruvate dehydrogenase by NADH and ACCOA) has important effects.
The overall control on integrated system behavior by NAD and ADP can be understood based on simulation of the model as follows. We define the rate of reducing equivalent (NADH and FADH 2 ) production as J DH ϭ J pdh ϩ J isod ϩ J akgd ϩ J sdh ϩ J mdh and compute the predicted steady-state J DH (6). For these simulations, the transport fluxes are omitted, and the tricarboxylic acid cycle is simulated in isolation. For this system, made up of reactions 5-15 in Table 1, the first 10 reactions sum to an overall reaction for the tricarboxylic acid (Reaction 1).
Since we do not account for transport into or out of mitochondria here, the reactants of the overall reaction of Equation 2 are clamped. In addition, aspartate and glutamate are held fixed, because there is no source or sink for these metabolites other than the glutamate oxaloacetate transaminase reaction. Since the electron transport system is not included, proton transport is not accounted for, and pH is held fixed. ) Since both ADP and NAD are not only regulators but also substrates for the tricarboxylic acid cycle, J DH vanishes when ADP and NAD levels go to zero. We can see from Fig. 8 that the relative NAD concentration is the more important controller of steady-state tricarboxylic acid cycle flux, and this point is consistent with the observation by LaNoue et al. (23). When ADP concentration is low, a variation in [NAD]/N o from 0 to 1 produces a significant change in normalized J DH from 0 to nearly 0.2. When NAD concentration is near zero, the rate of reducing equivalent production is not sensitive to ADP. However, the  flux is by no means insensitive to ADP; neither NAD nor ADP represents a sole independent controller of the system.
Regulation of ATP Synthesis by Inorganic Phosphate-Both experimental (25) and computational studies (5,32) show that the rate of mitochondrial ATP synthesis depends upon inorganic phosphate concentration. To investigate the role of phosphate as a regulator of cellular energetics, we examined the predicted phosphate metabolite concentrations in the integrated cell model while systematically removing mechanisms of energetic control by inorganic phosphate. First, we generated a model simulation of the in vivo relationship between work rate and phosphate metabolites under conditions identical to those used to obtain curve 1 in Fig. 7, with the exception that the P i concentration was clamped at the resting concentration (0.8 mM) for the purpose of computing the complex III flux. In other words, these model predictions, which are plotted as curve 2 in the figure, correspond to the case where phosphate activation of complex III is taken out of the model.
A third set of simulations of the integrated cell model was conducted with the matrix P i concentration clamped at 0.8 mM. The corresponding results are labeled curve 3 in the figure. Therefore, curve 3 represents predictions where all control related to changes in mitochondrial P i is removed from the model.
Both forms of the reduced model fail to reproduce the physiological response of the cell to changes in work rate. When regulation of complex III is not included, the predicted cellular ADP and P i concentrations are systematically higher than the measured data, and the energetic state of the cell (reflected in the ATP hydrolysis potential) is diminished. When matrix P i is clamped, the deviations from the experimental observations and the energetic state of the cell is more impaired compared with the normal case, because the effects of P i on both the tricarboxylic acid cycle activity and oxidative phosphorylation are not included in these simulations. Thus, inorganic phosphate concentration is a key signal in determining the mitochondrial response to cellular energy demands.
Effects of Cytosolic pH on Mitochondrial Function-Cytoplasmic pH in skeletal muscle tends to decrease during heavy exercise due to excess acidifying glycolytic flux (4,13). In addition, acidosis occurs in the heart during ischemia. To analyze the capacity of mitochondria to synthesize ATP during acidosis, the current model applied to in vivo skeletal muscle is simulated at different values of cytoplasmic pH. Plotted in Fig. 9 are model-predicted membrane potential and cytoplasmic ATP, ADP, and P i as functions of work rate for cytoplasmic pH values 6.4, 6.7, and 7.0. The predictions at pH 7.0 correspond to those reported in Fig. 7 for normal oxidative metabolism.
Decreasing the cytoplasmic pH results in a drop in mitochondrial membrane potential, which reduces the free energy level at which the ANT can deliver ATP to the cytoplasm. The result is a reduced concentration of ATP and increased concentrations of ADP and P i compared with normal. Thus, the model predicts that the oxidative work capacity of muscle decreases as the cytoplasmic pH value decreases.

DISCUSSION
Major Conclusions-The present work demonstrates that a vast amount of independent data, obtained from both in vivo and ex vivo systems, may be explained by a detailed model of mitochondrial energy metabolism. Predictions based on the parameterized and validated model suggest that the mitochondrial redox state is a primary regulator of tricarboxylic acid cycle flux. This model prediction is supported by a wealth of experimental observations (25,33,34). In addition, conclusions of previous studies that mitochondria redox state is strongly affected by available inorganic phosphate (5,25) are reinforced by the current study. The ability of our model to match in vitro data from preparations of isolated mitochondria and in vivo data from 31 P NMR spectroscopy in human subjects strongly depends on the influence of inorganic phosphate on tricarboxylic acid cycle kinetics and oxidative phosphorylation. Our analysis predicts that inorganic phosphate significantly influences tricarboxylic acid cycle flux through its roles both as a substrate and a necessary co-factor for transport of other substrates.
Based on the present modeling results and previous results focusing on oxidative phosphorylation, we propose that the control of mitochondrial ATP synthesis is dominated by ADPand P i -driven activation of oxidative phosphorylation and NAD-and P i -driven activation of the tricarboxylic acid cycle. Inorganic phosphate plays a significant role in stimulating both oxidative phosphorylation and tricarboxylic acid cycle, as originally proposed by Bose et al. (25).
Mitochondrial Metabolite Transporters-As illustrated in Fig. 1, the majority of tricarboxylic acid cycle intermediates are exchanged between the mitochondrion and its external environment. Thus, in order to use data from experiments using purified suspensions of mitochondria to identify a model of mitochondrial function that includes tricarboxylic acid cycle kinetics and to use the identified model to simulate in vivo function, it is necessary to account for the nine metabolite transporters and exchangers considered here. Since this model facilitates the simulation of mitochondrial suspensions respiring on different substrates under different conditions, it was possible to parameterize and challenge the model by a large number of independent experiments. According to the model parameterization, the majority of tricarboxylic acid cycle intermediate transporter fluxes operate near equilibrium. This is consistent with the observations of Williamson et al. (22,27) and justifies our use of simple mass action expressions for most transport fluxes. The AKG 2Ϫ / MAL 2Ϫ and ASP Ϫ /HGLU 0 exchangers are modeled using more detailed mechanisms to account for observed phenomena. The AKG 2Ϫ /MAL 2Ϫ exchanger flux (Equation B93 in Appendix B in the supplemental materials), which is modeled based on a rapid equilibrium random mechanism (35), competes with ␣-ketoglutarate, as observed by LaNoue et al. (26). The ASP Ϫ /HGLU 0 exchanger is unique in that it is driven by both the membrane potential and proton gradient (21,27,28). The ASP Ϫ /HGLU 0 flux expression (Equation B96 in Appendix B in the supplemental materials) is derived based on the rapid equilibrium random bi-bi mechanism with charge translocation, developed by Dierks et al. (21). Since both the AKG 2Ϫ / MAL 2Ϫ and ASP Ϫ /HGLU 0 exchangers are parts of the malateaspartate shuttle, their behavior is likely to be strongly coupled with substrate metabolism in the cytoplasm. In addition, several mitochondrial transport proteins have been shown to be nonselective and operate on multiple substrates. For example, the citrate carrier exchanges substrates citrate, isocitrate, and ␣-ketoglutarate, and the dicarboxylate carrier acts on succinate, malate, and inorganic phosphate (36,37). Because of similar structures and molecular weights among certain tricarboxylic acid intermediates, these intermediates may compete with each other for binding sites (19).
Thus, the mass action models used here for some of the metabolite transporters may be too simplified to mimic function under all physiologically relevant conditions. In particular, an oversimplified model for the tricarboxylate carrier (HCIT 2Ϫ /MAL 2Ϫ antiporter) is possibly responsible for the mismatch between simulations and experimental data in Fig.  2B. It is possible to improve the fit of Fig. 2B by lowering the activity of this transporter in the model. However, doing this increases the overall error for all experiments. Similarly, the simulated time scale of glutamate/aspartate exchange under state-2 conditions could be made to more closely match the experimental data in Fig. 2F by reducing the activity of the GLU Ϫ /H ϩ co-transporter. However, matching the data of Bose et al. (25), which were obtained using glutamate and malate as substrates, requires the relatively high GLU Ϫ /H ϩ co-transporter activity.
Challenges in Large Scale Computational Modeling-One major difficulty in constructing large scale integrated computational models of cellular biochemical systems is that they must be constructed based on components and data that are not always ideally compatible. For example, although the data used for identification and validation of the model presented in Figs. 2-6 of this study were all obtained from mitochondria isolated from heart, the developed model was then applied to explain data on skeletal muscle energetics in humans. Related to this issue is the fact that certain kinetic parameters used in this study were obtained from studies on enzymes obtained from different species and tissue types (e.g. kinetic parameters for citrate synthase were obtained from studies on enzyme obtained from rat liver and bovine heart; see section C3 in the supplemental materials). In addition, in this study, data of LaNoue et al. (23) obtained at 28°C were used in concert with data from Bose et al. (25) obtained at 37°C. As described in detail in the supplemental materials, the temperature effects on reaction thermodynamics were explicitly accounted for. However, since not enough data are available to develop a separate set of activity parameters for both data sets, the same enzyme activities are used for both temperatures. This fact may account for the systematic underprediction of the overall rate of oxidative phosphorylation measured at 37°C (see Fig. 4B).
Since issues of this sort are currently unavoidable in integrating models and data of the scale and scope addressed here, it is critical that such details are clearly and openly documented. As we have outlined under "Materials and Methods," we have exhaustively documented not only the sources of the kinetic parameter values used in this study in the supplemental materials, but we have also tabulated the species and tissue source along with alternative values where available.
Future Work-Although calcium has been shown to regulate mitochondrial energetics (38), the effect of calcium is not explicitly included in either oxidative phosphorylation or tricarboxylic acid cycle fluxes in this model. Our model prediction that substrate feedback primarily controls mitochondrial ATP synthesis in different energy states agrees with findings of Williamson, LaNoue, and other researchers (25,33,34). However, this finding does not conflict with the potential role of the calcium ion in allosteric regulation of certain tricarboxylic acid cycle enzymes (38). In fact, the role of calcium as a feed-forward signal in muscle cells has been demonstrated in previous modeling efforts (39). The parameterization of the current model is based on data where the calcium concentrations are expected to be saturating for allosteric binding calcium to regulatory elements and thus does not account for those calcium-regulatory mechanisms. Future applications of the current model will require that the role of Ca 2ϩ in regulating pyruvate dehydrogenase, isocitrate dehydrogenase, and ␣-ketoglutarate dehydrogenase be considered.
By providing a detailed description of mitochondrial metabolism that accounts for detailed biochemical thermodynamics, ion binding, and pH-dependent properties of biochemical reactions and ionic charge balance, the current model is the basis for future models accounting for additional pathways, including glycogenolysis and fatty acid metabolism. Such expanded models will be used to analyze experimental data sets and to investigate the regulation of energy metabolism in normal and diseased states. We propose that any extensions to the current model should be required to match the data used in the present study as well as or better than the current model does, in addition to matching additional data sets used to parameterize additional model components. Following this protocol will ensure that as computational models of cellular systems evolve in complexity and scale, their predicted behaviors are ideally matched to as much of the relevant available data as possible.