This Article Abstract Full Text (PDF) Appendix Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Cortassa, S. Articles by O'Rourke, B. Search for Related Content PubMed PubMed Citation Articles by Cortassa, S. Articles by O'Rourke, B.

An Integrated Model of Cardiac Mitochondrial Energy Metabolism and Calcium Dynamics

Sonia Cortassa^*, Miguel A. Aon^*, Eduardo Marbán^*,, Raimond L. Winslow^,* and Brian O'Rourke^*

The Johns Hopkins University, ^* Institute of Molecular Cardiobiology and The Center for Computational Medicine and Biology, Baltimore, Maryland

Correspondence: Address reprint requests to Brian O'Rourke, Ph.D., The Johns Hopkins University, Institute of Molecular Cardiobiology, 720 Rutland Ave., 844 Ross Bldg., Baltimore, MD 21205-2195. Tel.: 410-614-0034; Fax: 410-955-7953; e-mail: bor{at}jhmi.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
We present an integrated thermokinetic model describing control of cardiac mitochondrial bioenergetics. The model describes the tricarboxylic acid (TCA) cycle, oxidative phosphorylation, and mitochondrial Ca²⁺ handling. The kinetic component of the model includes effectors of the TCA cycle enzymes regulating production of NADH and FADH₂, which in turn are used by the electron transport chain to establish a proton motive force (µ_H), driving the F₁F₀-ATPase. In addition, mitochondrial matrix Ca²⁺, determined by Ca²⁺ uniporter and Na⁺/Ca²⁺ exchanger activities, regulates activity of the TCA cycle enzymes isocitrate dehydrogenase and -ketoglutarate dehydrogenase. The model is described by twelve ordinary differential equations for the time rate of change of mitochondrial membrane potential (_m), and matrix concentrations of Ca²⁺, NADH, ADP, and TCA cycle intermediates. The model is used to predict the response of mitochondria to changes in substrate delivery, metabolic inhibition, the rate of adenine nucleotide exchange, and Ca²⁺. The model is able to reproduce, qualitatively and semiquantitatively, experimental data concerning mitochondrial bioenergetics, Ca²⁺ dynamics, and respiratory control. Significant increases in oxygen consumption (V_O2), proton efflux, NADH, and ATP synthesis, in response to an increase in cytoplasmic Ca²⁺, are obtained when the Ca²⁺-sensitive dehydrogenases are the main rate-controlling steps of respiratory flux. These responses diminished when control is shifted downstream (e.g., the respiratory chain or adenine nucleotide translocator). The time-dependent behavior of the model, under conditions simulating an increase in workload, closely reproduces experimentally observed mitochondrial NADH dynamics in heart trabeculae subjected to changes in pacing frequency. The steady-state and time-dependent behavior of the model support the hypothesis that mitochondrial matrix Ca²⁺ plays an important role in matching energy supply with demand in cardiac myocytes.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
The high metabolic demand of the heart necessitates close matching of energy production to workload. Approximately 2% of cellular ATP is consumed on each heartbeat and almost all of this energy is replenished by mitochondrial oxidative phosphorylation under normoxic conditions (Harris and Das, 1991). Remarkably, the vital process of control and regulation of mitochondrial metabolism is not well understood.

After the development of the chemiosmotic theory of energy transduction by Mitchell (1961), several hypotheses for control of metabolism have been proposed. The first, referred to as respiratory control (Chance and Williams, 1956; Harris and Das, 1991), asserts that ADP availability to the ATP synthase is the limiting factor for mitochondrial ATP production. Increases in workload lead to a rise in cellular ADP and stimulation of oxidative phosphorylation, with the downstream "pull" on electron transport eventually causing acceleration of NADH production by the TCA cycle. The lack of a substantial increase in whole-heart ADP concentration in response to increased pacing frequency, and the observed decrease in µ_H produced by Ca²⁺ entry into mitochondria (Heineman and Balaban, 1990; Territo et al., 2000), have been used to argue against activation of ATP synthase based solely on its substrates (Harris and Das, 1991).

A second control hypothesis emerged from the discovery that Ca²⁺ activates mitochondrial matrix dehydrogenases (McCormack and Denton, 1984; McCormack et al., 1990). According to this hypothesis an increase in the workload of the heart is accompanied by a rise in cytosolic Ca²⁺, and the subsequent effect of matrix Ca²⁺ on the TCA cycle increases the supply of reducing equivalents (NADH, FADH₂) and an increased "push" of electrons through the respiratory chain, thereby increasing ATP production by generating more proton motive force. However, it has been argued that dehydrogenase activation by Ca²⁺ might not account quantitatively for observed changes in ATP synthesis rates (Harris and Das, 1991). Therefore, additional regulatory mechanism(s), including direct activation of the ATP synthase by Ca²⁺ (Territo et al., 2000), have been hypothesized to explain the acceleration of ATP production (Heineman and Balaban, 1990, Harris and Das, 1991). The present work addresses the hypothesis that Ca²⁺ is involved in coupling energy demand and supply in heart mitochondria. This hypothesis is tested through development and application of an integrative mathematical model that takes into account mitochondrial matrix- and membrane-based processes such as the TCA cycle, oxidative phosphorylation, and Ca²⁺ dynamics. To date, an integrated kinetic and thermodynamic model of cardiac energy metabolism incorporating mitochondrial Ca²⁺ dynamics has not been available. The present work builds a unified model of mitochondrial energy metabolism and examines its behavior with respect to changes in substrate supply, metabolic inhibition, and Ca²⁺. The results provide insight into the essential control properties of the system.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
General description
The model incorporates the main mitochondrial electrophysiological and metabolic processes and their interactions (Fig. 1). Model differential equations are given in Table 1, whereas rate equations and their behavior are presented in the Appendix. The major components of the model—mitochondrial matrix-associated processes, oxidative phosphorylation, and Ca²⁺ dynamics—are described in the following sections.

View larger version (35K): [in this window] [in a new window]   FIGURE 1  Schematic representation of the mitochondrial electrophysiological and metabolic processes and their interactions, as described by the model. The model takes into account oxidative phosphorylation and matrix-based processes in mitochondria. The tricarboxylic (TCA) cycle of the mitochondrial matrix is fed by acetyl CoA (AcCoA), which is the point of convergence for the oxidation of fatty acids and glucose, the two main substrates of the heart. The TCA (or Krebs) cycle completes the oxidation of AcCoA to CO2 and produces NADH and FADH2, which provide the driving force for oxidative phosphorylation. NADH and FADH2 are oxidized by the respiratory chain and the concomitant pumping of protons across the mitochondrial inner membrane establishes an electrochemical gradient, or proton motive force (µH), composed of an electrical gradient (m) and a proton gradient (pH). This proton motive force drives the phosphorylation of matrix ADP to ATP by the F1F0-ATPase (ATP synthase). The large m of the inner membrane (-150 to -200 mV; matrix negative with respect to the cytoplasm) also governs the electrogenic transport of ions, including the cotransport of ATP and ADP by the adenine nucleotide translocator, Ca2+ influx via the Ca2+ uniporter and Ca2+ efflux via the Na+/Ca2+ antiporter (Magnus and Keizer, 1997). The model also accounts for the explicit dependence of the TCA cycle enzymes isocitrate dehydrogenase and -ketoglutarate dehydrogenase on Ca2+. In this way, the rate of Ca2+ uptake by mitochondria is involved in membrane energization through the TCA cycle and oxidative phosphorylation, which are described in detail in the Appendix. Key to symbols: The concentric circles with an arrow across located at the inner mitochondrial membrane represent the m. Dotted arrows indicate regulatory interactions either positive (arrowhead) or negative (-).

The dashed lines in Fig. 1 show the regulatory feedback loops incorporated into the model. The rates of the highly-regulated dehydrogenases of the TCA cycle, namely, citrate synthase (CS), isocitrate dehydrogenase (IDH), -ketoglutarate dehydrogenase (KGDH) and malate dehydrogenase (MDH) and their regulatory interactions have been studied in detail and their behavior is presented in the Appendix. A fundamental feature of the present model is the explicit dependence of the dehydrogenases IDH and KGDH on mitochondrial calcium (Rutter and Denton, 1988; McCormack and Denton, 1984; Nichols et al., 1994; McCormack et al., 1990; Hansford and Zorov, 1998; also see Appendix).

Oxidative phosphorylation (OxPhos)
The OxPhos model is based on that of Magnus and Keizer (1997). The model describes NADH-driven electron transport and proton efflux, and F₁F₀-ATPase activity and its associated proton influx (Fig. 1; see also Appendix). The adenine nucleotide translocator and the proton leak are also taken into account. In general, this model follows a thermokinetic formulation. The driving forces are redox (A_res) and phosphorylation (A_F1) potentials (see Eqs. 26–34) and the proton motive force (µ_H). Coefficients in the equations correspond to rate constants of explicit conversion steps postulated in the Altman-King-Hill formulation of the model (Pietrobon and Caplan, 1985; Magnus and Keizer, 1997).

The F₁F₀-ATPase
The original expression for the rate of the F₁F₀-ATPase (Magnus and Keizer, 1997) was modified to account for reversibility of the enzyme when µ_H collapses (see Appendix). For example, in uncoupled mitochondria, hydrolysis rather than synthesis of ATP can occur and drive proton pumping.

Mitochondrial calcium dynamics
Mitochondrial Ca²⁺, [Ca²⁺]_m, is determined mainly by Ca²⁺ influx through the Ca²⁺ uniporter and the Na⁺-dependent Ca²⁺ efflux driven by the Na⁺/Ca²⁺ antiporter. The latter is assumed to be electrogenic based on the reported results of Baysal et al. (1994) (see Appendix).

Model building and numerical methods
The model was built using a modular approach. Briefly, each model component was examined separately, iteratively fine-tuning the kinetic parameters for the best fit to experimental data (Eqs. 13–39). Maple (Maple V, release 4, Waterloo, Canada) was used to study each analytical expression as a function of underlying parameters (see electronic Appendix: http://www.jhmi.edu/bor/model.htm).

Before the final model assembly, the TCA cycle and OxPhos modules were analyzed separately. The adequacy of the expressions for each of the processes involved was tested by simulating the model behavior at steady state. The system was considered to be in steady state when the magnitude of each time derivative was 10^-13. Model parameters such as the level of AcCoA, ion leak, and levels of electron carriers, were then varied to determine their influence on important mitochondrial functions. If the results were counterintuitive, or inconsistent with basic bioenergetic principles, the model was revised at the level of the analytical expressions stated in Maple and retested. The iteration process was performed until complete and satisfactory expressions were obtained, which converged to results consistent with available experimental evidence.

Numerical integration of model equations (Eqs. 1–12, Table 1) was performed using an Adams routine (SCoP, version 3.51w, Simulation Resources) until steady-state solutions (according to the criterion described above) were obtained. Steady-state values of each state variable were then used as input to software for performing bifurcation and continuation analysis (AUTO 1997, E. Doedel, Concordia University). This software was used to determine the dependence of steady-state solution properties (type and stability) on model parameters.

Dynamic simulations
To test whether or not the global kinetic response of the model accurately reproduces the dynamic response of mitochondrial energetics to changes in workload, we simulated the time-dependent behavior of mitochondrial NADH and Ca²⁺ during changes in the myocardial pacing frequency. Forcing pulse functions (SCoP version 3.51w, Simulation Resources) for the parameters [Ca²⁺]_i and [ADP]_i were incorporated into the model, and the simulations were compared with experimental data obtained in cardiac trabeculae (Brandes and Bers, 2002). Ca²⁺ transients were simulated with periodic square pulses of [Ca²⁺]_i to 1 µM (from a basal level of 100 nM) for 400 ms duration, applied at frequencies of 0.25 Hz or 2 Hz. To simulate the increased energy consumption when pacing frequency was increased, a single [ADP]_i pulse from 50 to 150 µM was applied for the whole interval of high frequency stimulation.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
Simulations of the response of mitochondrial energetics to changes in substrate input, levels of electron carriers, and cytoplasmic Ca²⁺ were performed to test whether the model qualitatively and semiquantitatively reproduced available experimental evidence. The overriding objective was to gain insight into the important sites and mechanisms of control of bioenergetics.

Control properties of the integrated model
Metabolic control analysis of oxidative phosphorylation has previously been used to demonstrate that control is distributed among the respiratory chain/substrate transport, the phosphorylation system, and the proton leak (Hafner et al., 1990; Murphy, 2001). This control can redistribute during electron chain inhibition, mitochondrial uncoupling, or with postischemic mitochondrial damage (Borutaite et al., 1995). Therefore, the simulations presented below examine the effects of metabolic perturbation on mitochondrial energetics to better understand control of the integrated system (Brand and Kesseler, 1995).

Influence of AcCoA on mitochondrial energetics
The contribution of substrate availability to the global response of bioenergetics was tested by varying either the earliest (AcCoA) or last (ADP) substrate (next section) of the integrated model. When the concentration of AcCoA was increased to a maximum of 10 µM, an approximate twofold increase in V_O2 was obtained; beyond this, V_O2 plateaued (Fig. 2 A). This increase in respiration rate is driven by the generation of NADH in the range of 0.001–1.5 mM (Fig. 2 A). The V_O2 registered at very low AcCoA is mainly driven by glutamate (see Appendix Eq. 24).

View larger version (20K): [in this window] [in a new window]   FIGURE 2  Effects of increasing substrate delivery on mitochondrial energetics. AcCoA concentration is varied and the steady-state behavior of the model analyzed according to a continuation analysis (see text). The correlative changes in the steady-state behavior of VO2 and NADH versus (A) AcCoA and (B) m are depicted when AcCoA is systematically changed. The simulation was run with the following set of parameters: VmaxANT = 5.0 mM s-1; pH = -0.6; res = 0.048 mM; ADPi = 0.05 mM; mM s-1; F1 = 0.06 mM; [Ca2+]i = 0.5 µM; gH = 0.01 mM s-1 V-1. The rest of the parameters are as indicated in Table 3. (C) Simulation results match experimental data obtained with pig heart mitochondria, redrawn from Territo et al. (2000).

Influence of phosphorylation potential on mitochondrial energetics
The respiratory rate (V_O2) increases as a function of the maximal rate of the adenine nucleotide translocator, ANT (Fig. 3 A) by six-to-sevenfold before saturation. In keeping with the classical theory of respiratory control, V_O2 varies with matrix ADP concentrations in the mM range (Fig. 3 A, inset). At low ADP, representing state 4 respiration, the basal rate of oxygen consumption is nonzero, supported in the model by a small proton leak across the mitochondrial inner membrane. A steep response of V_O2 to small changes in ADP, as would occur during increased metabolic demand, is evident as the transition is made from state 4 to state 3.

View larger version (15K): [in this window] [in a new window]   FIGURE 3  Relationship between membrane potential and respiratory flux upon changes in adenine nucleotide translocator activity. The changes in respiration rate, VO2, and membrane potential, m, are determined by the behavior of the whole model according to a steady-state analysis performed at different adenine nucleotide translocator maximal activities (VmaxANT) which result in different ADPm concentrations (A, inset). The left portion of curves a and b in panels A and B, corresponds to mitochondria in state 4 (i.e., low ADP) and the increases in VO2 result from the respiratory control exerted by increasing ADP concentrations. The simulations were run with the following parameter values: (a) pH = -0.6; res = 0.00072 mM; ADPi = 0.05 mM; mM s-1; F1 = 1.216 mM; [Ca2+]i = 0.1 µM; gH = 0.01 mM s-1 V-1; AcCoA = 1.0 mM. VmaxANT is varied in the range 0.01–25 mM s-1; (b) all parameters are the same as in a, except for res = 0.0051 mM and F1 = 0.61 mM. The remainder of the parameters are reported in Table 3. (C) Changes in the relationship between m and the respiratory rate obtained at various ratios of creatine/creatine phosphate in a preparation of heart mitochondria with normal or reduced respiratory activity (left-shifted curves). Redrawn from Borutaite et al. (1995).

View larger version (14K): [in this window] [in a new window]   FIGURE 5  Mitochondrial Ca2+ level after changes in cytoplasmic Ca2+. (A) The behavior of the whole model is studied with the AUTO 97 software as a function of cytoplasmic Ca2+ concentration, [Ca2+]i. The simulations were performed with the following set of parameters: (pull): RES = 0.00072 mM, F1 = 1.216 mM, s-1, s-1, and mM s-1, and (push): RES = 0.0036 mM, F1 = 0.0522 mM, s-1, s-1, and mM s-1. Curve a is obtained with the same parametric set as "pull" except for mM s-1. Otherwise the concentration of AcCoA is 1.0 mM, and the catalytic constants of CS and MDH are s-1 and s-1, respectively. (B) Experimental results redrawn from McCormack et al. (1990) showing the accumulation of mitochondrial Ca2+ in the absence (curve I) or in the presence of Na+ (curve II).

The results in Fig. 3 B are in qualitative agreement with the experimental results obtained with isolated heart mitochondria (reproduced in Fig. 3 C; Borutaite et al., 1995). Under conditions of a decrease in the respiratory activity of heart mitochondria due to ischemia, the curve describing the relation between _m and succinate-driven respiratory rate is shifted down and to the left with respect to the behavior of normal mitochondria (Fig. 3 C).

Influence of proton leak and respiratory inhibition on mitochondrial energetics
Inhibition of respiration is simulated by decreasing the concentration of electron carriers in the mitochondrial membrane (^res) (Fig. 4 A). The steepness of the V_O2 curve over the low range of ^res, and the electron carrier concentration at which V_O2 reaches saturation, are influenced by the magnitude of proton influx across the inner membrane, as depicted in Fig. 4 A, curves a and b. Proton influx can be varied by changing the concentration of F₁F₀-ATPase, which dissipates µ_H. Increasing the concentration of F₁F₀-ATPase increases V_O2 over the low ^res range (Fig. 4 A; curve a). As expected, with minimal ATP synthase activity, the membrane can support a higher _m at any V_O2 (Fig. 4 B; curve b). These modeling results resemble experimental data obtained with normal mitochondria in either state 3 or 4 upon titration with the metabolic inhibitor malonate at various creatine/creatine phosphate ratios (Fig. 4 C; Borutaite et al., 1995).

View larger version (13K): [in this window] [in a new window]   FIGURE 4  Relationship between mitochondrial inner membrane potential and the respiratory flux. The rate of the electron transport chain and the m are modulated through the concentration of electron carriers (res) (A, B). A steady-state analysis of the whole model behavior is performed. The parameters utilized in the simulations are those reported in Table 3 and as follows: Pi = 20.0 mM; pH = -0.6; mM s-1; ADPi = 0.05 mM; mM s-1; gH = 0.01 mM s-1 V-1, [Ca2+]i = 0.1 µM and AcCoA = 1.0 mM. The values of F1F0-ATPase concentration, F1, differ between curves, being 1.216 mM for curve a and 0.06 mM for curve b. (C) Changes in the relationship between m and the respiratory rate obtained upon malonate titration at different phosphorylation potentials (different ratios of creatine to creatine phosphate) in a preparation of heart mitochondria. Redrawn from Borutaite et al. (1995).

Steady-state behavior of mitochondrial calcium
The first step toward understanding the effects of Ca²⁺ on the global behavior of the model is to simulate the steady-state response of mitochondrial matrix free Ca²⁺ ([Ca²⁺]_m) to changes in free cytosolic Ca²⁺ concentration ([Ca²⁺]_i). Inasmuch as there is little quantitative information about the relative contributions of the two major Ca²⁺ transport pathways in determining steady-state [Ca²⁺]_m, the parameter space is explored with respect to the maximal Na⁺/Ca²⁺ antiporter rate. Results of these studies show that following [Ca²⁺]_i changes in the physiological range (0.01 to 1.2 µM), [Ca²⁺]_m increases nonlinearly after crossing a certain threshold [Ca²⁺]_i, and continues to increase to more than 20 µM without saturation (Fig. 5 A). This threshold depends on the ratio of Ca²⁺ uniporter and Na⁺/Ca²⁺ antiporter rates. In fact, the rightward shift of the [Ca²⁺]_m curve (Fig. 5, trace a) is due to a large increase in the rate of the Na⁺/Ca²⁺ antiporter, which counteracts influx through the Ca²⁺ uniporter. The general profile of [Ca²⁺]_m accumulation upon increasing [Ca²⁺]_i are in agreement with experimental results obtained by several authors (Fig. 5 B) (Miyata et al., 1991; Collins et al., 2001; reviewed by McCormack et al., 1990). The right shift in the experimental curve shown in Fig. 5 B was elicited by the addition of Na⁺ (McCormack et al., 1990) supporting the involvement of the Na⁺/Ca²⁺ antiporter in this response.

As the rates of the Ca²⁺ uniporter and Na⁺/Ca²⁺ antiporter are coupled indirectly through _m, the steady-state [Ca²⁺]_m uptake curve (Fig. 5 A) is relatively unaltered by changing other parameters that determine push or pull conditions. The effect of [Ca²⁺]_i on V_O2, however, is sensitive to such parameters, as described in the following section.

Response of mitochondrial energetics to changes in Ca²⁺
Effects on V_O2
Activation of V_O2 by Ca²⁺ is evident for both push and pull parameter sets (Fig. 6 A). Higher basal V_O2 values at a given [Ca²⁺]_i are evident for conditions in which the ANT and the electron transport flux are the main rate-controlling steps of respiration (pull) (Table 2). However, the degree of activation by Ca²⁺ is reduced compared to the push condition (9.1% versus 23%, Fig. 6 A) due to the Ca²⁺ insensitivity of the main rate-controlling steps of respiration. Paralleling the increase in V_O2, the steady-state fluxes through the Ca²⁺-sensitive enzymes, KGDH and IDH, also increases as a function of Ca²⁺ (Fig. 6 B, inset).

View larger version (15K): [in this window] [in a new window]   FIGURE 6  Respiratory flux and NADH levels upon changes in [Ca2+]I. The model is analyzed as described in Fig. 5 with the same set of parameters indicated in its legend for push and pull conditions. The steady-state respiratory flux (A) and NADH concentration (B) corresponding to the same simulations presented in Fig. 5 are depicted as a function of [Ca2+]i.

Effects on NADH and _m

View larger version (13K): [in this window] [in a new window]   FIGURE 7  F1F0 ATP synthase activity and membrane potential as a function of [Ca2+]i. The different curves represent the rate of the mitochondrial ATP synthase, VATPase (A) and m (B) at steady state after an increase in [Ca2+]i, as shown in Fig. 5 with the same set of parameters described in its legend (push and pull).

The net difference between the stimulatory effect of Ca²⁺ on TCA cycle dehydrogenases and the dissipative effect of Ca²⁺ on _m (Fig. 7 B) in turn determines the extent of stimulation or inhibition of ATP synthesis. Maximal stimulation occurs in the _m range of maximal responsiveness, i.e., near the midpoint of the V_ATPase versus _m curve (see Appendix Fig. A6, B). If _m is too high or too low, the increase in ATP synthesis rate is negligible, even if respiration is stimulated strongly. The same reasoning explains why, in a condition of pull V_ATPase is stimulated marginally (+5.2% from 24 nM to 1.07 µM [Ca²⁺]_i).

View larger version (14K): [in this window] [in a new window]   FIGURE A6  Dependence of F1F0-ATP synthase activity on ADPm, m and pH. The ATP synthetic (positive) or hydrolytic (negative) activity is plotted as a function of (A) ADPm concentration with pH = -0.4, and F1 = 0.525 mM, at different m as indicated (see Eq. 31); (B) m at the indicated pH values and a constant [ADP]m (=1.0 mM).

Under push conditions, the main rate-controlling steps of respiration are the Ca²⁺-sensitive dehydrogenases (KGDH and IDH) of the TCA cycle (Table 2). KGDH and IDH controlled 70% and 23%, respectively, of the respiratory flux when [Ca²⁺]_i stimulates the rates of respiration (Fig. 6 A), the F₁F₀-ATP synthase (Fig. 7 A) and the TCA cycle.

A large control by the respiratory chain and the ANT on the respiratory flux is obtained under pull conditions (Table 2). The latter situation corresponds closely to that described by Borutaite et al. (1995) regarding the control of respiration in heart mitochondria. Under these conditions, the processes downstream of NADH are flux-controlling (Table 2), but the ATP synthase is less stimulated by [Ca²⁺]_i (5.2%) when compared to the 18% increase registered under push conditions (Table 2). Thus, we conclude that [Ca²⁺]_i is better able to stimulate the rate of mitochondrial ATP synthesis when the TCA cycle exerts a significant control on respiration.

After a closer inspection of the initial reaction rate response of IDH and KGDH with respect to their substrates, only KGDH activity exhibits substantial elasticity to Ca²⁺ concentration within the stimulation range, in agreement with reported data (see Appendix Figs. A2 and A3; also see McCormack et al., 1990). This result explains why KGDH must exert a significant control on metabolism for respiratory flux to be stimulated by [Ca²⁺]_m (Figs. 5, and 6 A).

View larger version (14K): [in this window] [in a new window]   FIGURE A2  Isocitrate dehydrogenase activity as a function of substrate and Ca2+ concentrations. (A) The rate of IDH is studied as a function of isocitrate and Ca2+ concentrations in the physiological range (see Eq. 16). The analysis was performed with the following parameter values: [NAD+] = 9.5 mM; [ADP]m = 2.0 mM; mM; s-1 and those described in Table 3. (B) The IDH dependence on NAD+ concentration is analyzed at different ADP concentrations with [ISOC] = 0.5 mM, and [Ca2+] = 0.5 µM. Otherwise the parameters used are as in A.

View larger version (17K): [in this window] [in a new window]   FIGURE A3  Alpha-ketoglutarate dehydrogenase dependence on -ketoglutarate, NAD+, and Ca2+. (A) The set of parameters used in the analysis of the rate of KGDH (VKGDH) as a function of its substrate is as indicated in Table 3 (see also Eq. 17) with the following modifications: mM; s-1; [NAD+] = 1.0 mM; and [Mg2+] = 0.4 mM. (B) The analysis of the dependence of VKGDH on NAD+ concentration is performed with the following parameter values: [Mg2+] = 0.4 mM, and [Ca2+] = 0.1 µM. Otherwise the parameters used are as in A.

View larger version (32K): [in this window] [in a new window]   FIGURE 8  Time-dependent behavior of Ca2+ and NADH after changes in workload. (A, B) Experimental data showing the responses of NADH (A) and mitochondrial Ca2+ (B) to an increase in stimulation frequency from 0.25 to 2 Hz in a rat cardiac trabeculae (Reproduced from Brandes and Bers (2002) by copyright permission of the Biophysical Society). (C) Model response of mitochondrial Ca2+ to the pacing protocol described in Model Building and Numerical Methods, applied at low frequency (0.25 Hz) for 100 s, and then a high frequency (2 Hz) for 200 s, followed by a return to the low frequency. (D) Time course of NADH levels during changes in pacing frequency, simulated by increasing the Ca2+ pulse frequency combined with an increase in ADPi. The "pull" parametric condition is similar to that defined in the legend of Fig. 5 except for RES = 0.0006, F1 = 1.787 mM, and mM s-1; "pull II" is similar to "pull" except RES = 0.00072 mM and mM s-1. (E) The effects of a change in Ca2+ pulse frequency or ADPi applied separately.

Varying only the [Ca²⁺]_i transient driving function without changing ADP in the simulations results in a monotonic rise in NADH during an increase in pacing frequency (0.25 Hz–2 Hz) and a monotonic decline upon decreasing frequency (Fig. 8 E).

A pulse of ADP (50–150 µM) in the absence of an increase in [Ca²⁺]_i pacing frequency results in a rapid decrease in NADH due to stimulation of the respiratory chain (pull condition) and recovery upon returning to the basal level (Fig. 8 E).

When the change in Ca²⁺ transient frequency and the increase in ADP are combined, model NADH concentration transients closely resemble the experimental data (Fig. 8 D). The rate of recovery of NADH can be accelerated by a small (20%) increase in respiration and a fourfold increase in (Fig. 8 D; pull II). Fitting the curve of NADH recovery with a single exponential function renders t_1/2 values of 112 and 20 s for the conditions of pull and pull II, respectively. Mitochondrial Ca²⁺ relaxes with similar kinetics as NADH after the onset of high pacing frequency (t_1/2 = 127 and 33 s for pull and pull II, respectively), resembling more closely the relaxations observed by Brandes and Bers (2002).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
General
The main contribution of this modeling study is to link the TCA cycle, OxPhos, and mitochondrial Ca²⁺ dynamics through a mathematical model able to describe mitochondrial energetics based on well-known bioenergetic principles, i.e., the dependence of respiration on the redox potential and µ_H, and the F₁F₀-ATPase on the phosphorylation potential and _m. The model is used to explore the interaction between mitochondrial energy metabolism and important extramitochondrial signals such as Ca²⁺ and substrate delivery. Several features distinguish this model from those presented previously: 1), the model integrates the kinetics of the TCA cycle with oxidative phosphorylation and Ca²⁺ handling (Fig. 1), 2), both thermodynamics and kinetics of the reactions are taken into account, permitting dynamic studies of enzyme kinetics and key metabolic intermediates as well as redox potential and _m as the steady state is approached, and 3), the model incorporates the explicit dependence of the TCA cycle dehydrogenases on Ca²⁺, allowing examination of the influence of Ca²⁺ on mitochondrial metabolism as a whole.

Metabolic control of oxidative phosphorylation
The integrative framework provided by the model makes it possible to analyze points of control of mitochondrial respiration by investigating how mitochondrial ATP production is matched to an increase in cardiac workload through increases of cytoplasmic Ca²⁺.

Stimulation of respiration by ADP after changes in cellular ATP utilization has been proposed as one of the mechanisms by which mitochondrial ATP production meets the cellular ATP demand (respiratory control; see also Chance and Williams, 1956; Harris and Das, 1991). A shift from state 4 (respiration is low, µ_H is high, no ATP synthesis) to state 3 (respiration is high, µ_H is lowered by maximal ATP synthesis) is known to be accompanied by control redistribution in isolated mitochondria (Groen et al., 1982; Brown et al., 1990; Murphy, 2001). In cells, mitochondria are most often intermediate between states 3 and 4 (Ainscow and Brand, 1999; Murphy, 2001). The model reproduces the expected relationships between _m, V_O2, and F₁F₀ATPase rate after shifts between states 3 and 4 when experiments simulating inhibition of the respiratory chain or the adenine nucleotide translocator are performed (Figs. 3, 4).

Control redistribution also occurs during impairment of mitochondrial function after ischemic injury. Mitochondria isolated from ischemic hearts show a decrease in respiratory capacity and a decrease in phosphorylation potential when respiring in the presence of succinate. These effects are more pronounced the longer the ischemic period, e.g., from 30 to 45 min (Borutaite et al., 1995). Such an inhibition of the respiratory chain, simulated here by reducing the concentration of electron carriers, results in a decreased respiration rate and a concomitant decrease in membrane potential (Fig. 4), consistent with experimental results (Fig. 4 C; see also Borutaite et al., 1995).

The model also reveals a steep dropoff of oxygen consumption as the concentration of electron carriers (^res) reached a critical low level (Fig. 4 A). This type of threshold for disruption of mitochondrial function with respect to respiratory complex inhibition has been described previously (Murphy, 2001) and studied in vitro by titrating complex activity with irreversible inhibitors (Davey and Clark, 1996; Taylor et al., 1994; Rossignol et al., 1999). Maintenance of oxidative phosphorylation despite inhibition of a substantial fraction of the respiratory complexes may be explained by their high elasticities with respect to their substrates (Rossignol et al., 1999; Murphy, 2001). For example, partial inhibition of complex III will increase the ubiquinol/ubiquinone ratio, which would then increase the activity of uninhibited complexes. Thus, oxidative phosphorylation may be buffered against loss of respiratory chain components depending on the extent of inhibition of an individual complex required to decrease overall flux (Taylor et al., 1994). In our simulations, the threshold for collapse of respiration also depends on the inward flux of protons across the mitochondrial inner membrane (Fig. 4 A): when proton uptake is faster, the threshold ^res is shifted to a lower concentration range so that a greater fraction of carriers needs to be inhibited to affect respiration. It is important to note that ATP production would decrease before the decline in V_O2, as a result of the drop in _m at high proton uptake rates (Fig. 4 B).

Control by availability of substrate to the TCA cycle, AcCoA, may be relevant in the case of reperfusion of the heart after ischemia, inasmuch as TCA cycle flux is maintained or enhanced after brief ischemia (Weiss et al., 1993). The overall simulated flux through the TCA cycle in the range 0.05–0.1 mM s^-1 (10–20 µmol g dry wt^-1 min^-1), closely reproduces the global flux in whole hearts as reported in Jeffrey et al. (1999) and Williamson et al. (1976). Our results show that AcCoA levels exert control on respiratory flux, in parallel with NADH production, until saturation is reached (Fig. 2 A). As reported by Territo et al. (2001), substrate addition to isolated Ca²⁺-depleted mitochondria increases the respiration rate by 3.5-fold, _m by 30 mV (Fig. 2 C), and NADH by 7.7-fold. Similar trends are exhibited by V_O2, NADH, and _m in the simulations shown in Fig. 2.

Effects of Ca²⁺ on mitochondrial energetics
The inotropic mechanisms activated in response to an increase in cardiac workload (e.g., an increase in heart rate and/or Ca²⁺ release) also increase the time-averaged Ca²⁺ in the cytoplasm: this has been proposed as a control factor for increasing mitochondrial ATP production (Harris and Das, 1991; Territo et al., 2000, 2001).

An increase in [Ca²⁺]_i has two opposite effects on mitochondria: a stimulatory effect on the activity of the TCA cycle dehydrogenases as a consequence of higher [Ca²⁺]_m (Figs. 5 and 6), and a dissipative effect on the inner mitochondrial membrane potential due to the influx of the divalent cation. The balance between these two effects is reflected by _m. In isolated mitochondria, Ca²⁺ stimulation can increase NADH (McCormack et al., 1990; Hansford and Zorov, 1998; Territo et al., 2000), but the magnitude of this increase will also depend on any downstream effects on NADH oxidation. Flux through the electron transport chain will be accelerated by dissipation of µ_H, which could be caused by an increase in F₁F₀-ATPase velocity (due to increased ADP availability or direct stimulation of the ATPase (Territo et al., 2001)), or by non-ATPase-mediated energy dissipating pathways such as Ca²⁺ influx (relevant to the present model simulations) or enhanced proton leak. Non-ATPase energy dissipation will increase V_O2 without necessarily increasing ATP production; however, our model simulates both an increase in V_ATPase and an increase in V_O2 as a function of [Ca²⁺]_i (Figs. 6 A and 7 A), demonstrating that TCA cycle stimulation is adequate to bolster NADH levels and _m (Figs. 6 B and 7 B). The increase in the rate of ATP production by the mitochondria after an increase in [Ca²⁺]_i from 20 to 660 nM is 18%, a moderate change relative to the increase in V_O2 (usually more than 20% increase; Figs. 6 A and 7 A).

In the framework of the present model, the control of the respiratory flux at the mitochondrial level may be exerted by processes upstream and downstream of NADH (Table 2). Mitochondrial energetics is in a condition of push when respiration is controlled mainly by processes upstream of NADH, i.e., the TCA cycle. In this regard, a significant stimulation of ATP synthesis by Ca²⁺ can be achieved when respiration is controlled at the level of the Ca²⁺-sensitive dehydrogenases (Table 2: push). A drawback of the push condition is that the NADH levels produced by the model are lower than those determined experimentally (e.g., Brandes and Bers, 1997). Under pull conditions, downstream processes such as ANT or respiration itself, are the main rate-controlling steps of respiration (Table 2: pull).

Experimental data obtained for isolated rat liver cells show that control of respiration and oxidative phosphorylation is primarily downstream of NADH (Brown et al., 1990). Accordingly, rat liver mitochondria would be predominantly in the pull condition, inasmuch as 49% of the control is exerted by the processes of ATP synthesis, transport and consumption, 22% by proton cycling not coupled to ATP synthesis, and only 29% by respiration and upstream processes. Brown et al. (1990) estimated that 15–30% of the control over respiration would be exerted by processes involved in NADH generation. In our simulations, under pull conditions, the stimulation of oxidative phosphorylation by increasing cytosolic Ca²⁺ concentration is diminished (Figs. 6 and 7). Thus, the further incorporation of Ca²⁺-activated cytosolic ATPases, such as those involved in Ca²⁺ handling and muscle contraction, may be required in the future, to better couple the acceleration of V_O2 with the V_ATPase and increase the contribution of downstream pull on the system. In this regard, Territo et al. (2000) have reported that addition of Ca²⁺ to isolated mitochondria stimulated O₂ consumption as well as NADH levels. Determining ATP contents in mitochondrial preparations subjected to increasing Ca²⁺ concentrations, they inferred a stimulation of the ATP synthase activity. To account for these findings, they concluded that the F₁F₀-ATPase must be directly activated by Ca²⁺ (Territo et al., 2000, 2001). Nevertheless, it is notable that we can simulate an increase in ATPase activity upon increasing [Ca²⁺]_i without invoking other mechanisms of activation apart from the known biochemical and biophysical principles ruling mitochondrial bioenergetics. Stimulation depends only on the driving forces involved: the phosphorylation potential and the proton motive force. Obviously, inasmuch as the model simulations can reproduce different control properties of the system, they cannot provide an unequivocal quantitative answer to the relative contributions of upstream and downstream control mechanisms in living tissues. Therefore, we do not presume to strongly favor one control point over another for all cases; however, the model provides important feasibility arguments to be tested in future experimental studies.

An additional limitation of model interpretation lies in the absence of spatial organization of the mitochondria and cellular structures. For example, intracellular compartmentation of the adenine nucleotide pool has been proposed as an alternative explanation to Ca²⁺ activation in the context of matching of energy production and demand. It has been argued that cardiac contractile failure during ischemia or hypoxia still happens despite high ATP contents (70–80% of normal ATP content), leading to the hypothesis that a localized metabolic compartment may be present in which larger changes in phosphorylation potential occur. There is evidence for metabolic compartmentation (reviewed in Saks et al., 1998) and such metabolite channeling may be due to the close proximity of energy-consuming and energy-producing sites or to the actions of the creatine phosphate shuttle (Saks et al., 2001).

Similarly, the close apposition of sarcoplasmic reticular Ca²⁺ release sites and mitochondria may create a local Ca²⁺ microdomain, in which Ca²⁺ concentrations might be much higher than the bulk phase. Experimental evidence for fast uptake of mitochondrial Ca²⁺ has been reported (Rizzuto et al., 1998; Duchen et al., 1998). Our simulation results of the dynamics of Ca²⁺ uptake in response to pulses of [Ca²⁺]_i indicate that mitochondrial Ca²⁺ must essentially act as an integrator of cytosolic Ca²⁺ to reproduce intact muscle behavior (Fig. 8). This is another area in which spatial considerations and alternative Ca²⁺ uptake pathways may be incorporated in the future.

Time-dependent behavior of mitochondrial energetics
The ability of the model to simulate time-dependent behavior is indeed the most rigorous test of its validity, and emphasizes the utility of a model based on realistic kinetic properties of its components. We have been able to closely reproduce the NADH evolution in rat heart trabeculae when challenged by increases in workload through the pacing frequency of electrical stimulation (Brandes and Bers, 2002). The temporal profile of NADH evolution found experimentally can be simulated (compare Fig. 8, A and D). To reproduce the observed NADH transients, we have to incorporate an increase in ADP to represent increased cellular ATPase activity together with the increase in the pacing frequency of [Ca²⁺]_i. Thus, the model results suggest a co-participation of Ca²⁺ and ADP to explain mitochondrial NADH dynamics during an increase in workload. Additionally, we conclude that the under- and over-shoots of mitochondrial NADH during changes in pacing frequency are due to the rapid effects of ADP on respiration followed by the slower adaptation of the Ca²⁺-sensitive dehydrogenase activities. These results point out the involvement of respiratory control and activation of dehydrogenases by Ca²⁺ in the coupling between energy supply and demand in cardiac myocytes challenged by changes in workload.

Previous models of energetics
Previously, the most comprehensive metabolic models of cardiac metabolism were reported by Kohn and co-workers, to describe the energetics of the rat heart consuming different carbon substrates (Kohn and Garfinkel, 1979; 1983). Glycolysis, fatty acid oxidation, the TCA cycle, a lumped system representing mitochondrial respiration and oxidative phosphorylation, and chelation equilibria for nucleotides and organic acids were included in this model. However, these models were essentially based on mass action laws, with the assumption that most of the metabolic steps were reversible.

Other TCA cycle models were specifically developed for interpreting the data from NMR studies of whole hearts (Weiss et al., 1992, 1993; Yu et al., 1997; Tran-Dinh et al., 1996; Jeffrey et al., 1999). In such cases, the model equations were designed to fit the evolution of the ¹³C enrichment profiles of various TCA cycle intermediates, following the approach of Chance et al. (1983). The tracer equations used in these models do not address any explicit kinetic mechanism of each enzyme participating in the TCA cycle.

Other models of the TCA cycle have been developed for the amoeba Dictyostelium discoideum based on data about the temporal evolution of tracers (Wright et al., 1992; Shiraishi and Savageau, 1992). Wright and collaborators (1992) based their model on Michaelis-Menten kinetics taking into account the kinetic constants for substrates, products and some effectors obtained from the in vitro characterization of enzyme kinetics. Shiraishi and Savageau (1992), in turn, developed a model based on the power law formalism. Both models appropriately describe the behavior of the system for a certain range of parametric conditions. However, none of those models account for the effect of Ca²⁺, which is of central importance to cardiac mitochondrial energetics.

A model of pancreatic ß-cell mitochondrial metabolism by Magnus and Keizer (1997), which incorporates Ca²⁺ handling by mitochondria, involves three ordinary differential equations describing _m, Ca²⁺, and ADP concentrations as a function of time. This model takes into account the most important transport processes of the inner mitochondrial membrane, i.e., redox-driven proton pumps, the F₁F₀-ATPase, proton leak, the adenine nucleotide transporter, the Ca²⁺ uniporter, and the mitochondrial Na⁺/Ca²⁺ exchanger. Important limitations of this model are that it does not include stimulation of TCA cycle dehydrogenases by Ca²⁺ and that NADH is a fixed parameter (Magnus, 1995; Magnus and Keizer, 1997) instead of a variable (as in the present model). Nevertheless, in agreement with experimental data, the model behavior was able to reproduce the dependence of the rate of Ca²⁺ uptake by the uniporter on cytosolic Ca²⁺ concentration and predicts a sharp threshold for uptake at cytosolic Ca²⁺ concentrations of 0.4–0.5 µM (Magnus and Keizer, 1997).

An oxidative phosphorylation mathematical model has been developed by Korzeniewski (1996, 1998) and applied to the analysis of how the rate of ATP production by respiration is adjusted to meet energy demand during muscle contraction (Korzeniewski, 1998). As in previous cases, this model does not account for the coupling of OxPhos to the TCA cycle and thus no explicit dependence on Ca²⁺ by TCA cycle dehydrogenases was taken into account.

	CONCLUSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORMATION ACKNOWLEDGEMENTS REFERENCES

 
In summary, the integrated thermokinetic model constructed here reproduces the expected relationships among _m, V_O2, NADH, and mitochondrial matrix Ca²⁺ accumulation under various conditions including simulated inhibition of respiratory chain complexes and/or adenine nucleotide translocase.

Based solely on simple bioenergetic principles, it is demonstrated that Ca²⁺ activation of the TCA cycle, NADH production, mitochondrial respiration, and ATP synthesis is possible when respiratory flux is mainly controlled by the Ca²⁺-sensitive KGDH. During changes in workload, under pull conditions, both respiratory control and stimulation by Ca²⁺ appear to be required to reproduce the temporal evolution of NADH in heart trabeculae. Under other parametric conditions, supplemental control mechanisms for oxidative phosphorylation may be necessary to account for experimental data. Increased ATP production by Ca²⁺ can only be achieved when the extent of stimulation of NADH production exceeds the energy-dissipating effect of Ca²⁺ influx on _m.

This comprehensive model of cardiac mitochondrial metabolism will permit testing of bioenergetic control hypotheses, provide a basis for further refinement of individual steps of oxidative phosphorylation and mitochondrial Ca²⁺ handling, and allow for future incorporation of other biochemical pathways.

	APPENDIX

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX GLOSSARY SUPPLEMENTAL INFORM