HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 292/6/C2004    most recent 00271.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science 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 ISI Web of Science (4) Citing Articles via Google Scholar Google Scholar Articles by Nguyen, M.-H. T. Articles by Jafri, M. S. Search for Related Content PubMed PubMed Citation Articles by Nguyen, M.-H. T. Articles by Jafri, M. S.

CALL FOR PAPERS
Special Section On Mitochondrial Modeling and Function

Effect of Ca²⁺ on cardiac mitochondrial energy production is modulated by Na⁺ and H⁺ dynamics

¹Department of Bioinformatics and Computational Biology, George Mason University, Manassas, Virginia; ²Diagnostic Ultrasound, Woodinville, Washington; and ³Medical Biotechnology Center, University of Maryland Biotechnology Institute, Baltimore, Maryland

Submitted 17 May 2006 ; accepted in final form 1 March 2007

	ABSTRACT

TOP ABSTRACT THE MODEL NUMERICAL METHODS RESULTS DISCUSSION AND CONCLUSIONS GRANTS REFERENCES

 
The energy production of mitochondria in heart increases during exercise. Several works have suggested that calcium acts at multiple control points to activate net ATP production in what is termed "parallel activation". To study this, a computational model of mitochondrial energy metabolism in the heart has been developed that integrates the Dudycha-Jafri model for the tricarboxylic acid cycle with the Magnus-Keizer model for mitochondrial energy metabolism and calcium dynamics. The model improves upon the previous formulation by including an updated formulation for calcium dynamics, and new descriptions of sodium, hydrogen, phosphate, and ATP balance. To this end, it incorporates new formulations for the calcium uniporter, sodium-calcium exchange, sodium-hydrogen exchange, the F₁F₀-ATPase, and potassium-hydrogen exchange. The model simulates a wide range of experimental data, including steady-state and simulated pacing protocols. The model suggests that calcium is a potent activator of net ATP production and that as pacing increases energy production due to calcium goes up almost linearly. Furthermore, it suggests that during an extramitochondrial calcium transient, calcium entry and extrusion cause a transient depolarization that serve to increase NADH production by the tricarboxylic acid cycle and NADH consumption by the respiration driven proton pumps. The model suggests that activation of the F₁F₀-ATPase by calcium is essential to increase ATP production. In mitochondria very close to the release sites, the depolarization is more severe causing a temporary loss of ATP production. However, due to the short duration of the depolarization the net ATP production is also increased.

oxidative phosphorylation; tricarboxylic acid cycle; pH; computer model

Computational models have been valuable tools to help understand the complex dynamics of mitochondrial function in cardiac myocytes (3, 7, 20, 33, 42, 50). The contributions of these models will be discussed briefly. Aliev and co-workers developed a model of the compartmentalization of energy transfer from the mitochondria to the myofilaments that includes the compartmentalization of creatine kinase, and the kinetics of the creatine kinase reaction and ATP production (2, 3). The model demonstrates that energy transfer between functional compartments can play a role in the dynamics of cardiac energy metabolism. Beard (7) developed a model of the electron transport chain and the ATP synthase consisting of equations based on thermodynamic principles. The parameters of these equations were found by fitting the model outputs to experimental data. The model demonstrates that inorganic phosphate control of the electron transport chain is necessary to account for the experimental observations by Bose and co-workers (11). Cortassa and co-workers (20) integrated the Dudycha-Jafri model of the tricarboxylic acid (TCA) cycle with the Magnus-Keizer model for pancreatic -cell mitochondrial calcium dynamics. This model suggests that calcium can stimulate cardiac energy production. The Magnus-Keizer model includes the basic mechanism that regulated calcium as well as proton fluxes involved in maintenance of the membrane potential and ATP synthesis (46, 48). The model, however, is not based on the most recent data of cardiac calcium dynamics and does not balance protons or sodium. Jafri and Kotulska (33) developed a model of the cardiac mitochondria to explore the mechanisms of metabolic oscillations that depend on mitochondrial phosphate dynamics and inner membrane anion channels. This model successfully reproduces a large variety of experimental manipulations indicating the robustness of the formulation. Korzeniewski developed a model of mitochondrial energy metabolism based on empirical descriptions of various processes. However, it does not include the effects of calcium dynamics.

A computational model of mitochondrial energy metabolism has been developed to study how the various regulatory factors work together to control ATP production. In particular this manuscript explores how Ca²⁺ dynamics influences energy production by the mitochondria. The new model includes the following new features: 1) a new formulation of the calcium uniporter based on the experimental data of Kirichok and co-workers (41); 2) a new formulation of the Na⁺-Ca²⁺ exchanger (NCX) with a 3:1 Na⁺:Ca²⁺ stoichiometry as observed experimentally by Jung and co-workers (36); 3) a new formulation of the Na⁺-H⁺ exchanger (NHX) based on the experimental data of Kapus and co-workers (40); 4) a new formulation of the phosphate carrier based on the experimental data of Stappen and Kramer (56); 5) a new formulation of the F₁F₀-ATPase that includes a Ca²⁺ dependence based on the experimental work of Balaban and co-workers (4), a voltage dependence based on the experimental work of Kaim and Dimroth (5, 38), and a phosphorylation potential based on the experimental work of Westerhoff and colleagues (68); 6) proton buffering based on the experimental work of Vaughan-Jones and co-workers (63); and 7) description of regulation of the TCA cycle enzymes based on experimental kinetics (26). This will be discussed in more detail below.

As mentioned previously, this manuscript will focus on the contributions and mechanisms of Ca²⁺ regulation to cardiac energy metabolism. The model presented herein combines the Dudycha-Jafri model of the TCA cycle (26, 27) with the Magnus-Keizer model for calcium homeostasis, electron transport, and ATP production (46, 48). The Ca²⁺ transport across the mitochondrial membrane can be ascribed to two different mechanisms, Ca²⁺ uptake by the uniporter and Ca²⁺ extrusion by the NCX. The Magnus-Keizer model includes mitochondrial Ca²⁺ dynamics and its effects on energy metabolism in pancreatic -cells (46–48). It includes the following fluxes: Ca²⁺ uptake by the uniporter, Ca²⁺ extrusion by the NCX, the ADP/ATP translocase, a proton leak, ATP formation by the F₁F₀-ATPase, the respiration driven proton pumps, and H⁺ flux through the F₁F₀-ATPase. It contains differential equations describing the mitochondrial Ca²⁺ concentration ([Ca²⁺]_m), the mitochondrial ADP concentration ([ADP]_m), and the mitochondrial membrane potential (_m).

The Ca²⁺ uniporter in the Magnus-Keizer model is based on data from Wingrove and co-workers (69) from liver mitochondria, which shows steady-state mitochondrial Ca²⁺ dynamics. Furthermore, the model shows high steady state mitochondrial [Ca²⁺] (15.0 µM) in response to myoplasmic [Ca²⁺] (1.0 µM) that is uncharacteristic for cardiac myocytes. However, experimental works in cardiac myocytes by Wan and co-workers (65), Miyata and colleagues (49), and Chalmers and Nicholls (15) have shown more modest levels of steady state mitochondrial [Ca²⁺] in response to elevated myoplasmic [Ca²⁺].

Recent experimental developments have provided new data that can be used to improve the previous formulations for Ca²⁺ dynamics as well as provide the enzymes kinetics that are affected by Ca²⁺. Another recent study (13a) has suggested that mitochondrial Ca²⁺ dynamics can be quite fast, and refer to a second mode of Ca²⁺ transport named the "rapid uptake mode". In fact, Trollinger and co-workers (61) have observed beat-to-beat [Ca²⁺]_m transients that occur in conjunction with myoplasmic Ca²⁺ transients (see Figs. 11 and 12). They attribute the slow kinetics seen in earlier work to contamination of the fluorescence signal with measurement of lysosomal Ca²⁺ uptake and release. They observed that the calcium transients in the mitochondria are half the amplitude of the myoplasmic transients. Each has a similar rise time, but the decay is slightly slower in the mitochondria.

View larger version (32K): [in this window] [in a new window]   Fig. 1. A: schematic diagram of the mitochondria model. CoA, coenzyme A; FAD, flavin adenine dinucleotide; Pi, inorganic phosphate. B: schematic diagram of the tricarboxylic acid (TCA) cycle; m; number of protons needed to produce one ATP molecule.

View larger version (24K): [in this window] [in a new window]   Fig. 2. A: the Goldman-Hodgkin-Katz formulation was used to model the Ca2+ uniporter characterization by Kirichok and co-workers (41) ([Ca2+]e = 20 µM; [Ca2+]e = 100 µM; [Ca2+]e = 250 µM) using a voltage ramp from –160 to 20 mV (note they used the opposite sign convention) and mitochondrial Ca2+ concentration ([Ca2+]m) = 100.0 µM as in the experiments of extramitochondrial [Ca2+] ([Ca2+]e = 25 µM–long dashed line; [Ca2+]e = 100 µM–dashed line, and [Ca2+]e = 250 µM–solid line). In these simulations, the permeability of the Ca2+ uniporter (PCa) = 0.0006 s–1. m, Mitochondrial membrane potential. B: Na+ dependence of the Na+-H+ exchanger (NHX) based on experimental data () by Kapus and co-workers (40) and model fit (solid line). In these simulations, the maximal pump velocity (VNH) = 12 nmol/(mg/s) and the extra- and intramitochondrial pH is 6.95 similar to experiment. C: proton dependence of the NHX based on experimental data () by Kapus and co-workers (40) and model fit (solid line). In these simulations, VNH = 2.7 pH/s and the extramitochondrial pH is 7.0 similar to experiment. Different values for the maximal NHX flux were used to match the experimental results as in the experimental preparations there was likely different amounts of exchanger present. The units are different to match the experimental figures. D: F1F0-ATPase ATP production velocity based on the experimental data by Kaim and Dimroth (38) for Escherichia coli () and Propionigenium modestum (), and model fit (solid line).

View larger version (11K): [in this window] [in a new window]   Fig. 3. Simulated steady-state [Ca2+]m, resulting from different [Ca2+]e.

View larger version (22K): [in this window] [in a new window]   Fig. 4. A: simulated [Ca2+]m transients (solid line) in response to myoplasmic Ca2+ transients generated by the Jafri-Rice-Winslow (34) model of the guinea pig ventricular myocyte paced at 1.0 Hz (dotted line). B: simulated [Ca2+]m transients (solid line) at 2.0 Hz extramitochondrial Ca2+ transients (dotted line).

View larger version (18K): [in this window] [in a new window]   Fig. 5. Model responses to the pacing protocol from Fig. 4A. A: mitochondrial ATP concentration ([ATP]m: solid line) and mitochondrial proton concentration ([H+]m: dotted line). B: mitochondrial inorganic phosphate concentration ([Pi]m: dashed line), mitochondrial NADH concentration ([NADH]m: dotted line), and mitochondrial NAD+ concentration ([NAD]m: solid line). C: m (dotted line) and mitochondrial Na+ concentration ([Na+]m: solid line). The line style (dashed, dotted, or solid) is shown next to the vertical axis label for reference.

View larger version (19K): [in this window] [in a new window]   Fig. 6. Model ionic flux responses to the pacing protocol from Fig. 4A. A: Na+-H+ exchanger (NHX) flux (dashed line) and Ca2+ uniporter and Na+-Ca2+ exchanger (NCX) flux (overlapping traces: solid line and dotted line underneath, respectively). B: K+-H+ exchanger flux (long-dashed line), F1F0-ATPase proton flux (solid line) and respiration-driven proton pumps (dashed line).

View larger version (16K): [in this window] [in a new window]   Fig. 7. Mitochondrial energy metabolism flux responses to the pacing protocol from Fig. 4A. A: time derivative of the [NADH]m under normal (solid line) and "F1F0-ATPase Ca2+ dependence removal protocol" (dashed line). The two traces overlap. B: time derivative of [ATP]m under normal (solid line) and F1F0-ATPase Ca2+ dependence removal protocol (dashed line). C: m under normal (solid line) and F1F0-ATPase Ca2+ dependence removal protocol (dashed line). The two traces overlap.

View larger version (12K): [in this window] [in a new window]   Fig. 8. ATP production increases as a function of pacing frequency.

View larger version (21K): [in this window] [in a new window]   Fig. 9. Simulated Ca2+ marks at 1.0-Hz pacing frequency. A: Ca2+ transients near Ca2+ release sites (dashed line) and mitochondrial Ca2+ transients (marks) near the release sites (solid lines). B: m under normal (solid line) and voltage-clamp conditions (dotted line). C: time derivative of [NADH]m under normal (solid line) and voltage-clamp conditions (dashed line) during marks. The two traces overlap. D: time derivative of the [ATP]m under normal (solid line) and voltage-clamp conditions (dashed line) during marks.

View larger version (17K): [in this window] [in a new window]   Fig. 10. A: NADH production flux during 99.1% block of the Ca2+ uniporter (solid line) and during the combination of 99.1% block of the Ca2+ uniporter and 99.4% block of the NCX (dashed line) at various different extramitochondrial sodium concentrations ([Na+]e). B: [Ca2+]m during 99.1% block of the Ca2+ uniporter (solid line, left axis) at various different [Na+]m. The steady [Ca2+]m at different [Na+]e with no block (dot-dashed line, right axis) is much higher. The line style (dashed or solid) is shown next to the vertical axis label for reference.

View larger version (13K): [in this window] [in a new window]   Fig. 11. The effect of [Ca2+]e on steady-state mitochondrial energy production. A: ATP production flux by the F1F0-ATPase (Jp). B: NADH production flux by the TCA cycle (JNADHp). C: m.

Since the effectiveness of the NCX depends on the mitochondrial [Na⁺], the model must describe Na⁺ balance in the mitochondria. Sodium enters the mitochondrial matrix via the NCX. Sodium is extruded from the mitochondrial matrix across the inner mitochondrial membrane via the NHX, a passive pump that used the proton gradient to carry one proton into the mitochondria in exchange for one Na⁺ ion. The kinetics of the NHX has been characterized by Kapus and co-workers (39, 40), who demonstrated that the dependence of the NHX on H⁺ does not obey simple Michaelis-Menten type kinetics due to an allosteric binding site for mitochondrial H⁺. Kapus and co-workers explored the dependence of NHX on extramitochondrial Na⁺ as well as matrix pH in fluorescent BCECF-loaded rat heart mitochondria, a spectrally pH-sensitive fluorescent probe that allows for the simultaneous measurements.

The current mitochondrial models have not adequately addressed the issue of Na⁺ and H⁺ balance. For example, the Magnus-Keizer model assumed that although mitochondrial Na⁺ concentrations were not experimentally determined and since the NHX is fast, the pH gradient could inherently reflect the Na⁺ concentrations. To verify this, we have incorporated efflux and influx components to the model, that include the NHX, K⁺-H⁺ exchanger, phosphate carrier, and a pH buffer.

Another feature of the Magnus-Keizer model in need of modification is that a large Ca²⁺ influx would cause a fall in membrane potential that would shut down ATP production by the F₁F₀-ATPase, a process referred to as "short circuiting". In experiments, also accompanying the local myoplasmic Ca²⁺ transient, local depolarization of mitochondrial membrane potential has been observed, in agreement with the Magnus-Keizer prediction that elevated Ca²⁺ can depolarize the mitochondria (25). Others have attributed the "rapid uptake mode" to local elevated domains of Ca²⁺ and have measured mitochondrial [Ca²⁺] transients, named marks, in response to localized cytosolic Ca²⁺ transients called sparks (51). These local elevations of Ca²⁺ are likely to be much higher than the global elevations (16).

To address the new data on mitochondrial Ca²⁺ dynamics, this modeling effort has endeavored to incorporate the most recent available data on mitochondrial Ca²⁺ dynamics. To this end, the recent experimental characterization of the Ca²⁺ uniporter performed by Kirichok and co-workers (41) has been used to develop an accurate model of the uniporter. Accompanying this, a revised formulation of the NCX is also included, based on the experimental observations of Wingrove and Gunter, and modified to include the 3:1 stoichiometry of the exchanger found by Baysal and co-workers and Jung and co-workers (6, 36, 70, 71). Both of these formulations are used to replace the former descriptions to better simulate Ca²⁺ dynamics in the cardiac ventricular cell.

	THE MODEL

TOP ABSTRACT THE MODEL NUMERICAL METHODS RESULTS DISCUSSION AND CONCLUSIONS GRANTS REFERENCES

 
The mitochondrial model integrates the Magnus-Keizer model (46, 48) with a the model for the TCA cycle developed by Dudycha and Jafri (26, 27). It has also incorporated the following features: 1) an updated description of calcium dynamics to correspond to the cardiac ventricular myocyte; 2) proton balance; 3) sodium balance; 4) phosphate balance; 5) ATP balance; and 6) a new formulation of the F₁F₀-ATPase that describes experimentally observed dependence on calcium and membrane potential (Fig. 1A). Model parameter definitions and values are given in Tables 1, 2, and 3–5.

The TCA cycle model uses expressions based on Michaelis-Menten kinetics to describe the regulation of the key regulatory enzymes by various factors such as Ca²⁺, Mg²⁺, ADP, ATP, H⁺, NADH, and NAD⁺. The parameters in these expressions were simultaneously fit minimizing the sum of square error between the model equation and experimental data derived from multiple sources. With the use of an ANOVA F-test, parameters were tested to see whether they affected the fit in the model. Parameters that did not affect model results significantly were set to zero. Once the enzyme kinetics were defined, the enzymes were combined in a TCA cycle. Enzyme velocities were fit to match flux and steady-state concentrations based on experimental NMR glutamate enrichment data from Jeffrey and co-workers (35) and other biochemical assays in the literature. The details of the TCA cycle model follow.

The citrate synthase reaction converts oxaloacetate (OAA) and acetyl-coenzyme A (A-CoA) into citrate (Cit). The model (27, 28) fits to the experimental data indicates substrate inhibition by A-ACoA. The flux through citrate synthase (J_CS) can be described by

(1)

where K_M,A-CoA and K_S,A-CoA are binding constants for acetyl-CoA, K_M,OAA is the binding constant for oxaloacetate, K_i,A-CoA is the inhibition constant for acetyl-CoA, and V_CS is the reaction velocity.

Isocitrate dehydrogenase converts isocitrate (ISOC) and NAD⁺ into -ketoglutarate (KG) and NADH. Isocitrate dehydrogenase is activated by calcium and ADP and inhibited by NADH. The flux through isocitrate dehydrogenase (J_IDH) is described by

(2)

where K_M,ISOC is the binding constant for isocitrate, K_M,NAD is the binding constant for NAD⁺, K_H1, and K_H2 are binding constants for H⁺, K_a,ADP is the binding constant for ADP, K_Ca is the binding constant for Ca²⁺, and K_i,NADH is the inhibition constant for NADH. The exponents n_ISOC and n_NAD are the cooperativity factors for isocitrate and NAD⁺ binding.

The -ketoglutarate dehydrogenase complex converts -ketoglutarate (KG) and CoA into succinyl-CoA (S-CoA) and carbon dioxide producing NADH in the processes. It is activated by calcium and inhibited by S-CoA and NADH. The flux through -ketoglutarate dehydrogenase complex (J_KGDHC) is described by

(3)

where K_M,KG is the binding constant for KG, K_2,NAD is the binding constant for NAD⁺, K_d,Mg, and K_d,Ca are the dissociation constants for Mg²⁺ and Ca²⁺, V_KGDHC is the reaction velocity, and n_2,NAD is the cooperativity factor for NAD⁺ binding.

Malate dehydrogenase converts malate (MAL) and NAD⁺ into OAA and NADH. It is inhibited by oxaloacetate. The flux through malate dehydrogenase (J_MDH) is

(4)

where K_M,MAL is the binding constant for malate, K_i,OAA is the inhibition constant for oxaloacetate, K_3,NAD and K_4,NAD are binding constants for NAD⁺, K_H3, K_H4, K_H5, and K_H6 are H⁺ binding constants, k_offset is a fitting parameter that ensures that enzyme activity even at low pH, and V_MDH is the reaction velocity.

Succinate dehydrogenase converts succinate into fumarate. The reaction flux (J_SDH) is described by

(5)

where K_m,SUC is the binding constant for succinate, and K_i,FUM is the inhibition constant for fumarate.

The remaining reactions in the TCA cycle, aconitase (ACO), aspartate aminotransferase (AAT), and fumarate hydratase (FH), are described by mass action kinetics. For example, succinate-CoA ligase, which converts succinyl-CoA (SCoA) into succinate (SUC), can be described as

(6)

where k_f,SL is the forward rate constant and K_eq,SL is the reaction equilibrium constant.

This yields the following set of differential equations to describe the TCA cycle intermediates

(7)

where J_in is the flux of acetyl-CoA into the TCA cycle, J_ACO is the flux through aconitase, J_FH is the flux through fumarate hydratase, and J_AAT is the flux through AAT and the other fluxes are described above.

Mitochondrial Ca²⁺ dynamics is governed by three processes: 1) Ca²⁺ influx through the Ca²⁺ uniporter, 2) Ca²⁺ efflux through the NCX, and 3) Ca²⁺ buffering. One of the goals of the model is to simulate the recent experimental observations on mitochondrial Ca²⁺ dynamics in cardiac ventricular myocytes. A newly developed formulation of the uniporter describes the data observed by Kirichok and co-workers (41). The uniporter is assumed to be an ion channel permeable only to Ca²⁺. The model is based on a Goldman-Hodgkin-Katz formulation and closely matches membrane potential ramp data at different levels of myoplasmic [Ca²⁺] (Fig. 2A). The uniporter flux (J_UNI) can be described by

(8)

where P_Ca is the permeability of the Ca²⁺ uniporter, z_Ca is the charge of Ca²⁺, _m, and _e are the mitochondrial and extramitochondrial activity coefficients, and [Ca²⁺]_e is the extramitochondrial [Ca²⁺]. The fit of this formulation to the experimental data by Kirichok and co-workers is shown in Fig. 2A. In their experiment, they used pipette attached mitoplasts. The membrane potential was ramped from –160 to 20 mV (they used the opposite sign convention as Magnus-Keizer and us) and currents measured. The pipette solution contained no Ca²⁺ and EGTA to ensure that [Ca²⁺]_m remained zero. However, in reality, Ca²⁺ at the channel mouth is likely to be quite high. This is supported by the model, in which [Ca²⁺]_m was set to 100 µM to get a the proper reversal potential during the simulated ramp protocols.

The NCX in the Magnus-Keizer model shows a voltage dependence, which suggests that it generates current. However, it uses a ratio of two Na⁺ ions exchanged for one Ca²⁺ ion. Experimental results indicate that the stoichiometry is 3:1 for Na⁺:Ca²⁺ as observed by Baysal and co-workers and Jung and co-workers (6, 36). The NCX flux (J_NC) can be described by

(9)

where K_Na is the Na⁺ binding constant, K_Ca is the Ca²⁺ binding constant, and V_NC is the maximal exchanger velocity (33). This is based on a sequential mode of transport with random order of ion binding. The membrane translocation step is assumed to be voltage dependent as it is electrogenic. Rate constants for Na⁺ and Ca²⁺ binding are based on the experimentally measured values by Wingrove and Gunter (70), Crompton and colleagues (22), Coll and co-workers (17), and Paucek and co-workers (52).

The balance equations for [Ca²⁺]_m is described by

(10)

where _Ca is a constant describing the fraction of Ca²⁺ that binds to Ca²⁺ buffers in the mitochondria and t is time (in ms). This is approximated as 0.1 based on the data of Corkey and Williamson (19) and assumes that the binding of Ca²⁺ to buffers is fast compared with the other calcium fluxes.

The mitochondrial Na⁺ concentration ([Na⁺]_m) dynamics is governed by the NCX (described above) and NHX. Kapus and co-workers (39) observed that while the Na⁺ dependence of the Na⁺-H⁺ exchange obeyed Michaelis-Menten kinetics, the dependence on extramitochondrial proton concentration ([H⁺]_e) did not. They postulated the presence of a proton binding site on this exchanger outside the mitochondria. The following equation was derived to fit their data (Fig. 2, B and C) describing the Na⁺-H⁺ exchange flux (J_NH)

(11)

where K_H is the proton binding constant, K_Na is the Na⁺ binding constant, [Na⁺]_e is the extramitochondrial [Na⁺], and V_NH is the maximal pump velocity. The parameter K_H,reg, describing the binding constant for the allosteric site for mitochondrial H⁺, is empirically determined by fitting the data of Kapus and co-workers. The Na⁺ dependence of the exchanger is shown in Fig. 2B. The triangles are the experimental data of Kapus and co-workers (39) and the solid line the model fit. The proton dependence is shown in Fig. 2C, where the experimental data of Kapus and co-workers appears as squares and the model data appears as a solid line. The balance equation for [Na⁺]_m is described by

(12)

Proton (H⁺) balance in the mitochondria is governed by the F₁F₀-ATPase, the respiration driven proton pumps, the proton leak, the phosphate carrier, Na⁺-H⁺ exchange, K⁺-H⁺ exchange, and buffering. The formulation for the F₁F₀-ATPase (33) in the model that captures the dependence on the mitochondrial membrane potential (_m) based on the observation of Kaim and Dimroth (38), the phosphorylation potential described by Westerhoff and colleagues (68), and the dependence on [Ca²⁺]_m described by Balaban and co-workers (4). Experimental observations indicate that ATP production by F₁F₀-ATPase is mainly dependent on the mitochondrial membrane potential and not on the proton gradient (38). This is due to the properties that with large membrane potential, positive ions, such as protons, are driven into the mitochondria and can do so even up a concentration gradient. The formulation in the Magnus-Keizer model, which demonstrates a strong dependence on the proton gradient, was replaced by simpler formulation based on the work of Kaim and Dimroth (38). The model fit to their data is shown in Fig. 2D, where the squares show data from Escherichia coli and the triangles show data for Propionigenium modestum. The solid line shows the model fit. Similar dependence was found in cardiac mitochondrial by Berkich and co-workers (8). The rate of ATP production by the F₁F₀-ATPase (J_P) is described as

(13)

where K_ADP is the binding constant for ADP, K_ATP is the binding constant for ATP, K_Pi is the binding constant for P_i, K_Ca,ATP is the Ca²⁺ binding constant, and K_V,ATP is the membrane potential yielding half-maximal ATP production. The terms for the phosphorylation potential and Ca²⁺ dependence came directly from the experimental papers (4, 68). As there are three protons that pass through the F₁F₀-ATPase to produce one ATP molecule, the proton flux through the F₁F₀-ATPase (J_hf) is described as

(14)

The respiration driven proton pumps consume NADH to pump protons across the mitochondrial inner membrane out of the mitochondrial matrix. This process is also referred to as the electron transport or respiratory chain. The Magnus-Keizer model describes the electron transport chain as the single H⁺ efflux through the respiration-driven proton pumps (J_HR)

(15)

where the affinity of the driving reaction (A_res) is given by

(16)

and

(17)

This and other equations have been rescaled from nmol/min/(mg protein) to units of mM/ms using the following equation:

(18)

The rate of oxygen atom consumption (J_O) by the respiratory chain is described as

(19)

and the proton leak (J_HL) is

(20)

where g_H is the H⁺ leak conductance based on the Magus Keizer model.

The phosphate carrier carries phosphate into the mitochondria in one of three electroneutral transport modes: P_i/P_i counter-transport, P_i/OH^– countertransport, or P_i/H⁺ co-transport (56). The carrier is a bidirectional P_i/OH^– symporter that forms a ternary complex of P_i, OH^–, and the carrier. Therefore, the model is based on reversible Michaelis-Menten kinetics that involves two ternary complexes with a random order of addition of substrates resulting in simultaneous transport of the two substrates. The carrier (X) has the kinetic scheme where r_a, r_b, r_c1, r_c2, r₁, r₂ and r₃ are defined in Table 4.

in which P_i and OH can bind reversibly on either the cytoplasmic or mitochondrial side of the carrier in any order. When a P_i and an OH^– are bound on opposite sides of the carrier, a translocation can occur (middle). A mathematical description for this mechanism has been derived as follows

(21)

where V_PIC is the maximal transport rate, K_Pi,e is the extramitochondrial P_i dissociation constant, [P_i]_e is the extramitochondrial P_i concentration, K_OH is the OH^– binding constant, [OH^–]_m is the mitochondrial matrix OH^– concentration (calculated from [H⁺]_m), [OH^–]_e is the extramitochondrial OH^– concentration, [P_i]_m is the mitochondrial matrix P_i concentration, and K_Pi,m is the mitochondrial matrix P_i dissociation constant measured by Stappen and Krämer (56). When there is nonzero net P_i transport, i.e., the uptake and efflux rates are unequal, P_i/OH^– countertransport is most likely (56).

(22)

where [K⁺]_m is the mitochondrial [K⁺], [K⁺]_e is the extramitochondrial [K⁺], and X_KH incorporates exchanger velocity and a scaling factor. The balance equation for [H⁺]_m is described as

(23)

with _H being the proton buffering factor approximated as 0.00001 based on the work of Vaughan-Jones and colleagues (63) who estimated a reduction in the proton diffusion constant by a factor of 150,000 by buffering and the formula due to Wagner and Keizer (64) to convert free diffusion coefficients and buffering factor into an effective diffusion constant. The use of a constant buffering factor assumes that the buffers are fast compared with the other proton fluxes.

To capture the dynamics of mitochondrial phosphate ([P_i]_m) and [ADP]_m it is necessary to describe the phosphate carrier (above), the F₁F₀-ATPase (above), and the ADP/ATP translocase. The transport of ATP out of the mitochondria and ADP into the mitochondria is governed by the ADP/ATP translocase flux (J_ANT). The translocase exchanges one ATP for one ADP and is electrogenic. It is described by Magnus and Keizer as

(24)

where V_ANT is the velocity of the ADP/ATP translocase and f_p is the fraction of effective membrane potential for the translocase. The ionic forms of ADP and ATP are fractions of the ATP and ADP concentrations as defined by Magnus and Keizer (and Table 4).

The balance equations for [P_i]_m and [ADP]_m are

(25)

The remaining balance equations in the model describe mitochondrial NADH ([NADH]_m) and mitochondrial NAD⁺ ([NAD⁺]_m)

(26)

In the TCA cycle, NAD⁺ is converted to NADH by the enzymes isocitrate dehydrogenase (IDH), -ketoglutarate dehydrogenase complex (KGDHC), and malate dehydrogenase (MDH). These are consumed by oxidative metabolism in the electron transport chain. The final equation is for _m

(27)

In the model, the concentrations of extramitochondrial Na⁺, K⁺, H⁺, P_i, ADP, and ATP are assumed to be constant. This implicitly implies that the processes (binding and transport fluxes) that govern these concentrations are in equilibrium, i.e., they are instantaneously buffered.

	NUMERICAL METHODS

TOP ABSTRACT THE MODEL NUMERICAL METHODS RESULTS DISCUSSION AND CONCLUSIONS GRANTS REFERENCES

 
The 18 ordinary differential equations defining the model were solved using a fourth-order Runge-Kutta method with fixed time step of 0.02 ms. The computer code was written in FORTRAN and run on a Hewlett-Packard computer workstation (model C3000 running HP-UX or HP XW6200 running cygwin). Numerical results were visualized using PV-WAVE by Precision Visuals. Simulations involving steady state values were run for 60-s simulation time (by which time steady state was obtained). For simulations involving calcium transients, the calcium transients for both the cytosol and the subspace were generated by the Jafri-Rice-Winslow model (34). For simulation for pacing, the cytosolic calcium transients were used as the extramitochondrial calcium concentration. For simulations with marks, the subspace calcium concentration was used as the extramitochondrial calcium concentration. With pacing, the simulations were continued until fluxes and concentrations did not vary on a beat to beat basis.

Model parameters were found by three methods. They were determined finding experimental measured values in the scientific literature or calculations from these values. They were also found by fitting specific model equations to experimental data. The third method was to estimate the parameter so that model variable values and time courses matched experimental data. While this is in effect fitting the model to more macroscopic data this was done by hand rather than with a formal minimization algorithm. This third class includes rescaling parameters from previous models.

	RESULTS

TOP ABSTRACT THE MODEL NUMERICAL METHODS RESULTS DISCUSSION AND CONCLUSIONS GRANTS REFERENCES

 
The model was used to study both steady-state and transient changes to mitochondrial Ca²⁺ concentration ([Ca²⁺]_m). Parameter values from Tables 1 to 5 were used unless otherwise mentioned. First, mitochondrial variables were observed for conditions simulating a resting cardiac myocyte with extramitochondrial Ca²⁺ concentrations ([Ca²⁺]_e) set to 0.1 µM and extramitochondrial Na⁺ concentration ([Na⁺]_e) set to 5.0 mM. The steady-state values were close to experimentally observed quantities (Table 6). Reference values are also given with their source.

Figure 5 shows that the changes in myoplasmic [Ca²⁺] during pacing at 1.0 Hz lead to other changes in the mitochondrial model for energy metabolism. In Fig. 5A, [ADP]_m (solid line) declines with each [Ca²⁺] transient while [ATP]_m (dashed line) increases with each transient. The [H⁺]_m (dotted line) also oscillates in phase with the [ADP]_m oscillations, however, they are very small in amplitude (0.1 nM in amplitude) due to the high degree of proton buffering present in the mitochondria. The inorganic phosphate concentration ([P_i]_m) increases with each transient (Fig. 5B: dashed line) while the mitochondrial [NAD⁺] and [NADH] oscillate with [NAD⁺] (Fig. 5B: solid line), rising during the [Ca²⁺] transient when [NADH] (Fig. 5B: dotted line) falls (and vise versa) so that their sum remains constant. In Fig. 5C, _m (dotted line) depolarizes with each [Ca²⁺] transient. [Na⁺]_m increases during each [Ca²⁺] transient.

During the myoplasmic [Ca²⁺] transient, Ca²⁺ enters the mitochondria through the Ca²⁺ uniporter, producing a depolarizing current. This rise in [Ca²⁺] activates the NCX that extrudes Ca²⁺ ions and imports Na⁺ ions with a 3:1 ratio. This also produces a depolarizing current and is consistent with the depolarization of _m seen during the Ca²⁺ transient. Plotted in Fig. 6A is the Ca²⁺ uniporter flux (solid line) that overlaps the NCX flux (dotted line: not visible) and shows transient elevations during the extramitochondrial Ca²⁺ transients. Note that the Na⁺ entry is, in fact, three times the actual NCX Ca²⁺ flux shown and is similar to the Na⁺-H⁺ exchange flux (Fig. 6A: dashed line). During each Ca²⁺ transient [Na⁺]_m increases indicating that the Na⁺ extrusion through the NHX (Fig. 6A: dashed line) lags behind the Na⁺ entry through the NCX (Fig. 6A: dotted line but not visible). Accompanying each Ca²⁺ transient is increased activity of respiration-driven proton pumps, as indicated in Fig. 6B (dashed line) which serves to reduce [H⁺]_m. The other flux that decreases [H⁺]_m is K⁺-H⁺ exchange (long-dashed lines) which remains relatively constant. The flux of the F₁F₀-ATPase (solid line) increases [H⁺]_m and undergoes transient stimulation during the Ca²⁺ transients. This stimulation of the respiration driven proton pumps is primarily due to the transient decreases in _m (as is demonstrated in Fig. 9) as the decrease in _m lessens the electrochemical gradient that the proton pumps are working against. These transient decreases _m result from increased J_UNI and J_NC and to a lesser extent J_P during the Ca²⁺ transient. In summary, the transient decreases in membrane potential stimulate the respiration driven proton pumps which decreases [H⁺]_m. The decrease in [H⁺]_m causes increased NHX extrusion of Na⁺, which decreases [Na⁺]_m. The decreased [Na⁺]_m and large membrane potential allow efficient extrusion of Ca²⁺ by NCX.

The model can also be used to study how energy production by the mitochondrial is modulated during the Ca²⁺ transients. Figure 7A (solid line) shows the net NADH production flux (production-consumption). During each Ca²⁺ transient, there is increased consumption of NADH as indicated by the negative value of the trace. The trace returns to a positive value between Ca²⁺ transients. If the activation of the F₁F₀-ATPase by Ca²⁺ is removed by using the F₁F₀-ATPase Ca²⁺-dependence removal protocol protocol (setting [Ca²⁺]_m to the diastolic value in the F₁F₀-ATPase equation) there is little change to the net NADH production (dashed line). However, there is an impact on the ATP production rate (Fig. 7B). In the control case (solid line) there is a large transient increase in ATP production during each Ca²⁺ transient. However, in the simulation with F₁F₀-ATPase Ca²⁺-dependence removal protocol, there is no significant stimulation of ATP production. Figure 7C shows that the transient depolarization of _m is not significantly changed during the "F₁F₀-ATPase Ca²⁺ dependence removal protocol" (dashed lines) as it overlaps the control (solid). These simulations suggest that the Ca²⁺ transient stimulates energy production in parallel to the contraction induced by the Ca²⁺ transient.

The next set of simulations tests this hypothesis by exploring how mitochondrial energy production increases as pacing rate increases (Fig. 8). The Jafri-Rice-Winslow model was used to create myoplasmic Ca²⁺ transients at 1.0, 2.0, and 3.0 Hz. The average energy production of a 10-s period was determined (with [Na⁺]_e = 10.0 mM). For 0.0 Hz, a resting myoplasmic Ca²⁺ of 0.1 µM was applied. Once the model was at steady state, the average energy production rate was determined. As pictured, the energy production increases as pacing increases. One should note that this includes only the stimulation by Ca²⁺. It is likely that with increased workload, other factors such as changing NADH/NAD⁺ and ADP levels will also serve to modulate energy production. For example, if [ADP]_m is elevated by augmenting the rate of the ATP/ADP translocase, the degree of energy production potentiation increases with faster pacing to a greater extent.

Experimental studies have suggested that some mitochondria are also located near the Ca²⁺ release sites, i.e., near the cardiac dyadic subspace. In the case of Ca²⁺ sparks, localized mitochondrial [Ca²⁺] transients known as marks have been observed (51) accompanied by local depolarization of the mitochondrial membrane potential (25). The model was used to assess how these local Ca²⁺ signaling events affect energy metabolism. There is no experimental data measuring [Ca²⁺] in mitochondrial near the dyadic subspace during the cardiac cycle. This being the case, these studies are predictions of the model. To do so, the mitochondrial model was exposed to the elevated Ca²⁺ likely to be seen near the dyadic subspace generated from the Jafri-Rice-Winslow model for the guinea pig ventricular myocyte. Figure 9A shows the elevated local myoplasmic [Ca²⁺] near the subspace (dashed line) and the resulting mitochondrial [Ca²⁺] (solid line). The resulting transient decrease in _m in Fig. 9B (solid line) shows large depolarization similar to the local depolarization observed by Duchen and co-workers (25). To demonstrate the role of _m during the Ca²⁺ transients, _m was held constant in a voltage clamp protocol (dotted trace). Figure 9C (solid line) shows that a transient increase of net NADH production (production by TCA cycle minus consumption by the respiration driven proton pumps) followed by decrease in net NADH production occurs during depolarization. If there is no depolarization (dotted line), there is a large transient increase in the net NADH production due to activation of the Ca²⁺-dependent dehydrogenases in the TCA cycle (it is no longer masked by the increase in NADH production due to depolarization). This is followed by a large transient decrease below the baseline as the proton pumps catch up and exceed the production of NADH. The transient increase in NADH production by the TCA cycle is also seen in the control case (solid line) by the notch. Figure 9D shows that during the [Ca²⁺] transients, there is also a transient increase in the ATP production flux (solid line). However, during the control protocol, there is a marked reduction in ATP production in the early transient due to the depolarization of the membrane potential (solid line). When the membrane potential is clamped, the transient inhibition is relieved (dotted trace). Despite this transient inhibition of ATP synthesis, the average ATP production is still augmented by 29% in these mitochondria over the mitochondria located away from the subspace.

To verify that the model correctly predicts the dependence of energy metabolism on mitochondrial Ca²⁺ dynamics, the experiments of Cox and Matlib (21) were simulated. They performed experiments to assess the role of the NCX on oxidative energy production in isolated rabbit heart mitochondria. They measured NADH production in the presence of 1.0 µM ruthenium red to block Ca²⁺ uptake by the uniporter and in the presence of 1.0 µM ruthenium red and 10.0 µM CGP-37157 to block both Ca²⁺ uptake by the uniporter and Na⁺-Ca²⁺ exchange. They observed NADH production at various [Na⁺]_e ranging from 0.0–10.0 mM. With ruthenium red alone, as [Na⁺]_e increased from 0.0 to 5.0 mM, NADH production rate dropped to one-half its value at 0.0 mM Na⁺. When both ruthenium red and CGP-37157 were added, the NADH production rate declined modestly as the extramitochondrial [Na⁺] increased from 0.0 to 10.0 mM. The model simulates these experiments as shown in Fig. 10A. When uptake by the Ca²⁺ uniporter was decreased by 99.1%, the NADH production rate falls with increasing extramitochondrial [Na⁺] in a similar fashion as the experiments. When the uniporter was decreased by 99.1% and Na⁺-Ca²⁺ exchange 99.6% blocked, the NADH production rate remained declined similar to experiment. The mechanism behind this reduction is shown in Fig. 10B. As extramitochondrial [Na⁺] concentration increases the [Ca²⁺]_m is declines. This decline with increasing sodium (solid line, left axis) is qualitatively similar to the decline observed by Wan and co-workers (65). For comparison, the dependence of [Ca²⁺]_m in the case of no block is shown (dot-dashed line, right axis). It is much higher than the case with block. The decline of [Ca²⁺]_m decreases activation of the Ca²⁺-sensitive dehydrogenases (IDH and KGDHC) in the TCA cycle which decreases NADH production. It also decreases ATP production by the F₁F₀-ATPase, which results in less depolarization of the membrane potential, which leads to less NADH consumption. This allows NADH to increase with further inhibits the TCA cycle.

Figure 3 shows that steady-state mitochondrial [Ca²⁺] increases as extramitochondrial [Ca²⁺] rises similar to experiment. Figure 11 explores the impact this has on energy production. In Fig. 11A, the ATP production flux initially rises to peak at 0.2 mM [Ca²⁺]_e and then declines. While this seems contrary to the idea that Ca²⁺ activates energy metabolism it is consistent with the oxygen consumption results by Wan and co-workers (65) who see similar biphasic behavior with a peak at 0.2 µM [Ca²⁺]_m. The model suggests that the mechanism behind this. The TCA cycle NADH production flux (Fig. 11B) shows similar biphasic behavior; however, it plateaus at a much higher level indicating that is not the reason for the decline. The membrane potential, however, declines montonically (Fig. 11C). As the membrane potential is crucial for ATP production by the F₁F₀-ATPase its decline results in a decline in ATP production. The biphasic behavior results because at low Ca²⁺ activation of the TCA cycle and F₁F₀-ATPase by Ca²⁺ can more than compensate for the decline in membrane potential. The decline in membrane potential results from increased depolarizing current through the Ca²⁺ uniporter and Na⁺-Ca²⁺ exchange.

	DISCUSSION AND CONCLUSIONS

TOP ABSTRACT THE MODEL NUMERICAL METHODS RESULTS DISCUSSION AND CONCLUSIONS GRANTS REFERENCES

 
A model for mitochondrial energy metabolism in the cardiac ventricular myocyte has been developed. It improves upon the Magnus-Keizer mitochondria by addition of cardiac Ca²⁺ handling mechanism, Na⁺ balance, proton balance, phosphate balance, ATP balance as well as model of the TCA cycle. There have been differences in calcium handling between cardiac ventricular and pancreatic -cell mitochondria observed experimentally (61, 65). Furthermore, there might be some other functional differences due to the functional differences between mitochondria in the two cell types. For example, mitochondrial in cardiac ventricular myocytes are very active in producing ATP, while the metabolic demands of the pancreatic -cell are less. Furthermore, mitochondria in the pancreatic -cell seem to play a role in glucose sensing for insulin production (45, 53).

The model produces Ca²⁺ dynamical responses consistent with experimental data under both steady-state Ca²⁺ levels and simulated Ca²⁺ transients during pacing. The steady-state Ca²⁺ data closely matches the results of Wan and co-workers (65) both quantitatively and quantitatively. The transient data closely simulates the experimental results of Trollinger and co-workers (61). Similar to their experiments, increasing extramitochondrial Na⁺ decreases mitochondrial Ca²⁺. It is interesting to note that in both the experimental and simulated data steady-state mitochondrial [Ca²⁺] exceed extramitochondrial [Ca²⁺] unless the concentration is small (i.e., 0.1 µM). During calcium transients the mitochondrial [Ca²⁺] is always smaller than the extramitochondrial [Ca²⁺]. This occurs because the transients are brief so that the mitochondrial [Ca²⁺] does not reach its steady-state value given the rate of Ca²⁺ uptake by the uniporter.

The work also explores the contributions of Ca²⁺ dynamics to the regulation of energy production and the importance of parallel activation of the TCA cycle and the F₁F₀-ATPase. Increases in extramitochondrial Ca²⁺ lead to increases in mitochondrial calcium. This activates energy metabolism both at the TCA cycle through the isocitrate dehydrogenase and -ketoglutarate dehydrogenase complex and at the F₁F₀-ATPase. If the F₁F₀-ATPase is not activated by Ca²⁺, ATP production cannot be stimulated by Ca²⁺. This occurs because accelerated TCA cycle flux leads to an increase in NADH concentration (and a decrease in NAD⁺), which feeds back to inhibit the TCA cycle at enzymes such as isocitrate dehydrogenase. Acceleration of the F₁F₀-ATPase, Ca²⁺ uniporter, and NCX depolarizes the mitochondrial membrane resulting in increased activity of the respiration driven proton pumps and increased consumption of NADH.

In addition to the simultaneous activation of ATP production, the term parallel activation has been used by Korzeniewski (42) to describe the observation that both ATP usage and production were both stimulated by calcium. ATP production would be increased by increased TCA cycle activity producing more NADH and hence more of a proton gradient. Increased contraction leads to increased ATP utilization, which would increase ADP, which could then also activate mitochondrial energy production. The current model does not include ATP utilization as that is left for future work.

The model shows an increase in energy production from a simulated resting state up to a simulated 3.0-Hz pacing protocol. The energy production increases by a factor of 1.84 going from rest to pacing at 3.0 Hz. While this is a significant increase, there are other factors that might contribute to increased energy production as the heart rate increases. Experiments by Brandes and Bers (12, 13) show a transient drop in cellular NADH autofluorescence when the pacing rate is increased in a stepwise fashion. They see the opposite with decreased pacing. Because most of the cellular NADH is in the mitochondria, this indicates that mitochondrial NADH likely drops resulting stimulation of the TCA cycle. Changes in ADP might also serve to regulate mitochondrial energy production.

While the model incorporates advancements over previous models there are additional improvements possible. For example, in the model the electron transport chain is represented by a single equation. This can better be represented by modeling the activity of each cytochrome individually as was done in the model by Beard (7). Incorporation of such a model will better represent the metabolism of NADH and FADH₂ as FADH₂ consumption is not easily accounted for in the present formulation. Furthermore, it will allow easier model comparison to the enzyme-dependent fluorescence recovery after photobleaching experiments which can measure enzyme activity of NADH production (18). Another possible improvement is the inclusion of other ion channels such as inner membrane ion channels or the ATP-sensitive K⁺ channel. Finally, the formulation of the K⁺-H⁺ exchanger while taken from a published work is not e