Regulation of mitochondrial respiration in heart cells analyzed by reaction-diffusion model of energy transfer

Marko Vendelin1, Olav Kongas1, and Valdur Saks2,3

1 Institute of Cybernetics and 2 Laboratory of Bioenergetics, Institute of Chemical and Biological Physics, Tallinn, Estonia; and 3 Joseph Fourier University, BP 53X-38 041 Grenoble Cedex, France


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

The purpose of this study is to investigate theoretically which intracellular factors may be important for regulation of mitochondrial respiration in working heart cells in vivo. We have developed a model that describes quantitatively the published experimental data on dependence of the rate of oxygen consumption and metabolic state of working isolated perfused rat heart on workload over its physiological range (Williamson JR, Ford G, Illingworth J, Safer B. Circ Res 38, Suppl I, I39-I51, 1976). Analysis of this model shows that for phosphocreatine, creatine, and ATP the equilibrium assumption is an acceptable approximation with respect to their diffusion in the intracellular bulk water phase. However, the ADP concentration changes in the contraction cycle in a nonequilibrium workload-dependent manner, showing the existence of the intracellular concentration gradients. The model shows that workload-dependent alteration of ADP concentration in the compartmentalized creatine kinase system may be taken, together with the changes in Pi concentration, to be among the major components of the metabolic feedback signal for regulation of respiration in muscle cells.

compartmentation; adenosine diphosphate; creatine kinase; metabolic oscillations; mathematical modeling


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

THE CLASSICAL PARADIGM of cellular energy metabolism, i.e., its conceptual framework, is based on the assumption of the equilibrium (or at least quasi-equilibrium) of the reactions involved. Practically all experimental studies of muscle metabolism have been performed during the last two decades within the framework of these theoretical considerations, with creatine kinase (CK) equilibrium used for calculation of cytoplasmic ADP concentration and all derived thermodynamic parameters (23, 34). However, these theoretical concepts are based on very few experimental works exclusively performed on resting muscles, where the equilibrium state is the easily expected one (34). Recent reinvestigation of this problem with a mathematical model of compartmentalized energy transfer has shown that this concept is not totally valid: in the working heart cells the CK reaction is clearly out of equilibrium during most of the contraction cycle (1, 32). This observation raises the following question: For which compounds involved in the energy metabolism of the cell are the calculations based on the assumption of equilibrium reasonably correct? Also, another equally important question relates to the existence and deepness of the concentration gradients of ADP and other metabolites in the cells (20). Answers to these two questions are important for understanding the nature of the intracellular factors regulating the rate of mitochondrial oxidative phosphorylation in the cells in vivo. We investigated these problems theoretically by analysis of different versions of mathematical models of energy transfer in muscle cells developed on the basis of the previous model of Aliev and Saks (1). This model was supplemented with the model of oxidative phosphorylation developed by Korzeniewski and Froncisz (21, 22) to account for the changes in mitochondrial membrane potential and to calculate the rates of oxygen consumption (VO2) for different workloads. The new model describes quantitatively the experimental dependencies of metabolic parameters of isolated perfused working rat hearts on the workload over the physiological range of VO2 and workloads published by Williamson et al. (38). The results of analysis of intracellular metabolic changes confirm that there are significant oscillations of the cytoplasmic ADP concentration in the cells within the cardiac cycle and significant concentration gradients of this metabolite in the working muscle cells because of the nonequilibrium state of the CK reaction. It is concluded that localized changes of ADP concentrations (its oscillations) due to the nonequilibrium mode of functioning of the CK system in the cells open new possibilities of the metabolic regulation of mitochondrial respiration, in comparison with the equilibrium state of the cellular metabolism.


    GLOSSARY
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

General


AK Adenylate kinase
ANT Adenine nucleotide translocase
CK Creatine kinase
Cr Creatine
IM Intermembrane
PCr Phosphocreatine

Model variables


Met Concentration of metabolite (ATP, ADP, AMP, PCr, Cr, Pi) in myofibril and cytoplasm
Meti Concentration of metabolite in mitochondrial IM space
ATPg, ADPg ATP and ADP concentrations in microcompartment
UQ Oxidized form of coenzyme Q
c3+ Oxidized form of cytochrome c
NAD+ NAD+ in mitochondrial matrix
Hx H+ in mitochondrial matrix
ATPx ATP in mitochondrial matrix
Pi,x Pi in mitochondrial matrix

Diffusion


 nabla 2 Laplace operator
DMet Diffusion coefficient of metabolite
GCK Spatial distribution function of myofibrillar CK
GAK Spatial distribution function of myofibrillar AK
GH Spatial distribution function of myofibrillar ATPase

Myofibrillar ATPase


HATP Rate of ATP hydrolysis in myofibril
HATP(max) Maximal HATP

CK in myofibril, myoplasm, and IM space


vCK CK reaction rate
vMiCK vCK in IM space
vMiCK,G vCK, ATP and ADP from microcompartment
vMiCK,I vCK, ATP and ADP from IM space
V1, V-1 Maximal CK reaction rates in forward and reverse directions
KiaKbKicKdKidKibKIb Dissociation constants for MgATP, Cr, MgADP, PCr, PCr, Cr, and Cr from CK-MgATP, CK-MgATP-Cr, CK-MgADP, CK-MgADP-PCr, CK-PCr, CK-Cr, and CK-MgADP-Cr complexes
KiaG, KicG MiCK dissociation constants for ATP and ADP in microcompartment

AK in myofibrils and IM space


vAK AK reaction rate
vAKmit vAK in IM space
kfakba AK forward and backward reaction rate constants

Mg2+ buffering


mMet Mg2+-bound forms of metabolite (ATP or ADP)
fMet Mg2+-free forms of metabolite (ATP or ADP)
Mgext Mg2+ in cytoplasm and myofibril
Mgx Mg2+ in mitochondrial matrix
KDText, KDDext, kDTx, kDDx Mg2+ dissociation constants

H+ buffering


Hext, Hx H+ concentration in cytoplasm and matrix
pHext, pHx pH in cytoplasm and matrix
 Delta pDelta pH, Delta psi Protonmotive force, pH difference, and membrane potential across inner mitochondrial membrane
rbuff, rbuff,0 H+ buffering capacity coefficients

Diffusion through mitochondrial outer membrane


<A><AC>n</AC><AC>&cjs1164;</AC></A> Unit vector, outward normal to the domain boundary
Li Width of the layer between inner and outer membranes
Lm Width of the layer between core of myofibril and mitochondrial outer membrane
RMet Permeability of the mitochondiral outer membrane for the metabolite

ATP diffusion between microcompartment and IM space


vexchATP, vexchADP Rate of ATP or ADP flux between microcompartment and IM space
RexchATP, RexchADP Rate constants

ATP/ADP translocation by ANT


vANT Net rate of ATP export by ANT from matrix
vright-arrow  GANT Rate of ATP export by ANT from matrix to microcompartment
vright-arrow  IANT Rate of ATP export by ANT from matrix to IM space
VANT Maximal net rate of ATP export by ANT from matrix
KgKi ANT reaction dissociation constants

Fractional volumes


FVi Fractional volume of IM space with respect to cell volume
FVx Fractional volume of matrix with respect to cell volume

Substrate dehydrogenation


vdh, kdh Net and maximum rates of substrate dehydrogenation
KmN Michaelis-Menten constant of substrate dehydrogenation for the NAD+/NADH ratio
pD Relative sensitivity coefficient of substrate dehydrogenation to the NAD+/NADH ratio

Phosphate carrier


vPI, vf,PI, vb,PI Net, forward (to matrix), and backward rates of phosphate carrier
pKa One-half dissociation constant for monovalent Pi

Proton leak


kL1, kL2 Phenomenological constants for proton leak

ATP synthesis


vsn, ksn Net and maximal rates of ATP synthase
 Delta G0, Delta Gsn Gibbs free energy and thermodynamic span of ATP synthase
na H+/ATP stoichiometry for ATP synthase

Respiratory chain


vC1, vC3, vO2 Rates of complexes I, III, and IV
kC1, kC3, kC4 Maximum rates of complexes I, III, and IV
a2+, c2+, a3+, c3+ Concentrations of reduced and oxidized forms of cytochromes a3 and c
kMO Michealis-Menten constant of complex IV for oxygen
En, Eu, Ec, Ea Redox potentials of NAD, ubiquinone, and cytochromes c and a3
En,0, Eu,0, Ec,0, Ea,0 Standard redox potentials of NAD, ubiquinone, and cytochromes c and a3
Ax, NADtot, O2 UQtotctotatot Total concentrations of adenine nucleotides, NAD, oxygen in matrix, ubiquinone, and cytochromes c and a3.
rbuff, rbuff,0 Buffering capacity coefficients

General constants


T, R, F Absolute temperature, gas constant, and Faraday number
Z, u, eta Constants for unit conversion


    METHODS OF MODELING
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

Principles. The spatially inhomogeneous reaction-diffusion model of energy transfer considers the reactions in three main compartments of cardiac cells: the myofibril together with the myoplasm, the mitochondrial intermembrane (IM) space, and the mitochondrial inner membrane-matrix space (Fig. 1). The metabolites described by the model in the myofibrils and IM space are ATP, ADP, AMP, phosphocreatine (PCr), creatine (Cr), and Pi. All these metabolites diffuse between the cytosolic and IM compartments, where the metabolites are involved in the CK and adenylate kinase (AK) reactions. In addition, the ATP is hydrolyzed in the myofibrils. In the IM space the mitochondrial CK reaction is coupled to the adenine nucleotide translocase (ANT) reaction via strong coupling between these enzymes; the coupling is moderated by a diffusional leak of the intermediates. The metabolites described by the model in the matrix compartment and in the inner membrane are NADH, coenzyme Q, cytochrome c, protons, ATP, ADP, and Pi. Three coupled reactions representing the production of protonmotive force by complexes I, III, and IV are included in the model. Protonmotive force is consumed by ATP synthase and membrane leak. The ANT rate is considered to depend on membrane potential. Pi is transported by a phosphate carrier.


View larger version (30K):
[in this window]
[in a new window]
 
Fig. 1.   Scheme of 2-dimensional compartmentalized energy transfer model of cardiac cells. A: components of model. Total diffusion path in transverse direction is 1.2 µm consisting of mitochondrial outer membrane (heavy line below mitochondrion), 0.2 µm for myoplasm, and 1.0 µm for myofibril. Bottom edge of A corresponds to core of myofibril. Total diffusion path in longitudinal direction is 1.0 µm. Creatine kinase (CK) and myosin magnesium ATPase activities are distributed nonuniformly along myofibril (37). B: activities of enzymes in dependence of longitudinal coordinate. CK concentration in myoplasm is taken to be equal to that in I-band.

Model description. The mathematical modeling was started by reproducing the one-dimensional model (1) of compartmentalized energy transfer from mitochondria to the center of myofibrils in the cardiac cells. It was then converted to a spatially inhomogeneous two-dimensional reaction-diffusion model (Fig. 1). ATP utilization by myofibrillar ATPases and regeneration from PCr and ADP by CK and AK in the myofibril and the myoplasm are considered two-dimensional processes. Namely, the active ATPases are not distributed uniformly along the myofibril but are found only in the regions where actin and myosin overlap (Fig. 1). According to Wegmann et al. (37), the myofibrillar CK isozyme is also distributed inhomogeneously within the myofibril. The distribution functions for CK and AK in myofibrillar space were constructed on the basis of Fig. 4 in Ref. 37 and are shown in Fig. 1B. The distributions of CK and AK in the myoplasm are uniform, and their activities are assumed to be equal to those in the I-band. Then the two-dimensional model was supplemented with oxidative phosphorylation processes by combinination with the model by Korzeniewski (21).

The new model was used to calculate the concentration of ADP and its gradients in the cell within the contraction cycle, concentration of metabolites, and VO2. During contraction, the myofibrils are supposed to hydrolyze ATP, with kinetics predicted from change in the time derivative of pressure in isovolumic rat hearts, i.e., linear increase in ATP hydrolysis rate up to 30 ms followed by its linear decrease to zero at 60 ms at a heart rate of 333 beats/min (180 ms for the cardiac cycle). The total amount of ATP hydrolyzed by myofibrillar ATPases is used to regulate the workload from zero to its maximum value of up to 960 µmol ATP · g dry wt-1 · min-1 corresponding to ~160 µmol O2 · g dry wt-1 · min-1 for working heart (32, 37). We use also a spatially homogeneously distributed basal level of ATP consumption in the myofibrillar and myoplasmic compartment with 18 µmol ATP · g dry wt-1 · min-1 corresponding to ~3 µmol O2/ · g dry wt-1 · min-1. Addition of basal level consumption is of some importance at very low workloads. The lengths of the full diffusion paths are taken to be 1.0 µm in the "longitudinal" direction and 1.2 µm in the "transverse" direction to the myofibril (5). These are the distances between the Z and M lines and between the myofibril core and the mitochondrial outer membrane, respectively. The radius of the myofibril is taken to be 1.0 µm. Thus in the model the myoplasm covers the area 1.0 × 0.2 µm (Fig. 1). The fractional volumes of the cell compartments for myofibrils, myoplasm, IM space, and matrix are chosen to be 10:2:1:3. We assume that within the mitochondrial IM space and matrix the concentrations of the metabolites depend only on the longitudinal coordinate. In the model of Korzeniewski (21), we modified slightly the constants in the membrane leak function to meet the following two conditions: 1) at 40% of maximal VO2, 25% of oxygen is consumed by proton leak (Fig. 3 in Ref. 14 and Fig. 1 in Ref. 10) and 2) for arrested rat heart, the minimal VO2 is 12-15 µmol O2 · g dry wt-1 · min-1 (32).

To study the importance and significance of the use of the two-dimensional model of compartmentalized energy transfer, we compare it with one- and zero-dimensional models. The dimension of the model is changed by variation of the diffusion coefficients. To ignore a longitudinal or transverse dimension, the corresponding diffusion constant is increased by 105 times. Thus the diffusion in this direction is considered to be infinitely fast and equalizes the concentration of the metabolites, thus resulting in zero partial derivatives along this direction. In the one-dimensional case, the longitudinal (along myofibrils) dimension is ignored and the myofibrillar part of the model, except for AK, is reduced into the original model by Aliev and Saks (1). In the zero-dimensional case, both dimensions are ignored. To study the reduced diffusion in the myofibrillar compartment (3), we reduced the diffusion coefficients by a factor of 10. To investigate the equilibrium case for metabolites, we increased the activity of CK 105-fold compared with its normal activity. It is supposed that in this case the reaction rates establish the equilibrium fast enough compared with the diffusion rate. In these different models we also study the influence of the outer mitochondrial membrane on the transport processes (see APPENDIX). Functional coupling of mitochondrial CK (miCK) to ANT was assumed to occur by means of high local ATP and ADP concentrations in a 10-nm (1, 13) narrow-space microcompartment between coupled enzymes. We have made the following assumptions in our description of the ANT and miCK. First, ANT is translocating the adenine nucleotides between matrix space and microcompartment and, partly, IM space. Second, the miCK adenine nucleotide binding center may be occupied by the adenine nucleotides from the microcompartment or from the IM space. Furthermore, it is assumed that when the reaction involves ATP from the microcompartment, the miCK reaction product ADP will be released also into the microcompartment. This kind of channeling of ADP from miCK to ANT ("reverse coupling") was not considered in the previous model (1). Third, diffusion between the microcompartment and the IM space is restricted. Finally, the capacity of the microcompartment is considered infinitely small; therefore, for each adenine nucleotide the influx is taken to be equal to its efflux. With these assumptions, the local ADP and ATP concentrations are determined by the ANT reaction, the miCK reaction, and the restricted diffusion. The values of the coefficients of the diffusion restriction from the microcompartment and the dissociation constants of adenine nucleotides from miCK were estimated from 1) deviation of the mass action ratio of the miCK from the equilibrium constant value observed experimentally (31, 32) in coupled rat heart mitochondria under conditions of oxidative phosphorylation and 2) change in the apparent dissociation constants of miCK due to the oxidative phosphorylation (19).

All equations included in the model are shown in the APPENDIX. All values of the experimentally determined parameters used in the model were taken from literature and are shown in Table 1.

                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Values of experimentally determined parameters

Computation. The model equations were numerically solved by a finite-element method in conjunction with Galerkin's method. The resulting system of ordinary differential equations was solved by the backward differentiation formula that is able to treat stiff equations. The accuracy of the solution was tested by comparing different spatial discretizations and varying the tolerance of the ordinary differential equation solver. The finite-element discretization was performed using the software package Diffpack (7), and the system was integrated using the DVODE package (6). The calculation time varied from 1 s to a couple of hours on an Alpha-Linux PC, depending on the complexity of the model.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

Testing the new model for "in vitro" conditions. We first tested the new model for its ability to describe the published experimental data in vitro. Figure 2 shows that the new model describes adequately the experimentally observed increase in the rate of PCr production in the coupled miCK reaction under the condition of oxidative phosphorylation in isolated rat heart mitochondria described by Jacobus and Saks (19) and Saks et al. (31). The increase in the rate of PCr production by oxidative phosphorylation is explained by increased turnover of adenine nucleotides in the miCK and oxidative phosphorylation reactions coupled by ANT (31, 32). Another important observation is that, under conditions of oxidative phosphorylation, direct supply of ATP and removal of ADP by ANT exert strong kinetic and thermodynamic control of the miCK, making the reaction of PCr synthesis favorable (19, 31, 32). Under these conditions, inhibition of the oxidative phosphorylation results in rapid reversal of the miCK reaction, and production of ATP, rather than PCr, becomes favorable because of kinetic and thermodynamic properties of the isolated CK reaction (19, 31, 32). This alteration in the behavior of the miCK in the presence and absence of oxidative phosphorylation is quantitatively described by our new version of the model (Fig. 3). Thus the inclusion of Korzeniewski's model of mitochondrial reactions (21) and the detailed description of coupling between miCK and ANT used in the new version of the model give a satisfactory description of the mitochondrial reactions of energy production and their control in vitro.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 2.   Dependence of rate of phosphocreatine (PCr) production in mitochondrial CK reaction on ATP concentration in medium with () and without () oxidative phosphorylation (ox phos). Experimental points for isolated rat heart mitochondria are from Jacobus and Saks (19). Solid lines, calculations by using model described in this work. Reaction medium contained 5 mM Pi, 15 mM creatine (Cr), and 50 µM ADP.



View larger version (12K):
[in this window]
[in a new window]
 
Fig. 3.   Analysis of CK reaction in isolated rat heart mitochondria. Mitochondrial CK reaction was reversed by inhibition of oxidative phosphorylation with oligomycin. Conditions of reaction correspond to those used in experiments by Saks et al. (31): 2 mM PCr, 0.12 mM ATP, 0.05 mM ADP, 40 mM Cr, 5 mM Pi. Calculated curve for relative rate [rate of reaction related to maximal velocity (Vmax) of forward reaction of PCr synthesis in mitochondria] is shown. At time 0, oxidative phosphorylation was taken to be inhibited, and this quickly reversed reaction from PCr synthesis [positive CK reaction rate (vCK)] to its utilization (negative vCK), as shown in experiments by Saks et al.

Theoretical description of experimental results on isolated perfused hearts "in vivo." It is important to know exactly how the model describes the cardiac energy metabolism in vivo for the wide variety of conditions, notably for the workload changes. There is abundant information in the literature on the metabolic changes in the perfused heart during workload transitions, but in too many of these studies the workload is changed in a rather narrow range (usually 1.5- to 3-fold), so one can only guess what may really happen outside these limits (15, 17, 33). Therefore, these works are not suitable for serious modeling. Fortunately for us, there is one classical study by Williamson et al. (38) in which Neely's working isolated rat heart perfusion protocol was used to change the workload and, correspondingly, VO2 and energy fluxes >10-fold, up to probably maximal possible values. Also, changes in metabolite levels were very carefully analyzed in this work (38). For this reason, we used their experimental data as a basis of application of our model for analysis of the in vivo situation in the heart cells.

Figure 4 shows the calculated (our results) and experimental values [results of Williamson et al. (38)] of VO2 as a function of the workload that was altered by changing the rate of ventricular filling of isolated rat heart perfused according to the "working heart" protocol of Neely (38). The experimental and theoretical dependencies are nearly linear functions perfectly fitting each other. The maximal VO2 recorded by Williamson et al. (38) is close to 160 µmol O2 · g dry wt-1 · min-1 (38), which, according to the data of Mootha et al. (25), corresponds closely to the maximal activity of mitochondria calculated from state 3 respiration rates in vitro and tissue contents of mitochondria. The model describes equally satisfactorily the stable levels of main metabolites: ATP level always stays at its initial value, and PCr level starts to decline significantly only when VO2 exceeds 80-100 µmol O2 · g dry wt-1 · min-1. PCr-to-ATP and PCr-to-Cr ratios (very often the only parameters measured in an NMR experiment) change from their initial values of 1.5 very slowly at moderate workloads and more significantly at the highest workload (Fig. 5), in good accordance with the experimental data of Willamson et al. These results are in accordance also with all published data of Balaban et al. (2).


View larger version (15K):
[in this window]
[in a new window]
 
Fig. 4.   Computed (solid line) and experimental () O2 consumption rates (VO2) of working cardiac muscle. Experimental points, obtained in experiments with isolated rat hearts perfused according to working heart protocol of Neely, in which workload was changed by increasing ventricular filling rate, are from Fig. 5B of Ref. 38. Relative workload is fraction of maximal workload applied [maximal filling rate (38)]; in computations this is rate of ATP hydrolysis by actomyosin ATPase. Zero workload simulates situation in arrested heart, where O2 is consumed to sustain only basal level of ATP consumption and membrane leak through mitochondrial inner membrane. CK+, normal rates in hearts with intact CK activity. open circle  (VO2 = 41 µmol O2 · g dry wt-1 · min-1), computed maximal VO2 by CK-deficient hearts (CK-). Computed relation for CK-deficient hearts coincides with similar relation for normal heart at low workload. Relation is not absolutely linear, i.e., slope of curve at left edge is ~10% lower than that at right edge; reason for this nonlinearity is exponential dependence of membrane leak on protonmotive force. Calculations were made using 0-dimensional model.



View larger version (13K):
[in this window]
[in a new window]
 
Fig. 5.   Average PCR-to-Cr and PCr-to-ATP ratios over cardiac cycles as functions of VO2 calculated with 0-dimensional model. Model solutions are compared with experimental values of ratios of PCr to ATP () and PCr to Cr () from literature (32).

Under conditions of this remarkable metabolic stability, however, the model shows a very significant change in ADP during the contraction cycle in the myofibrillar space, and the peak values of ADP increase with elevation of workload and, correspondingly, of VO2 (Fig. 6). Minimal ADP levels are characteristic for the diastolic phase of the contraction cycle and close to the equilibrium values of ADP concentration. In accordance with many earlier observations, they do not change with elevation of workload and stay close to 50 µM (dotted line in Fig. 7). The average values of ADP concentration per cycle change more significantly (solid line in Fig. 7). At the same time, ADP peak value increases >40-fold when VO2 changes 10-fold (Fig. 7). To accommodate the traditional way of thinking in physiological and biochemical audiences, which are more accustomed to Michaelis-Menten-type analyses of reactions, in Fig. 8 these relationships are presented as respiration rate vs. ADP concentration. The relationships between VO2 and ADP concentration resemble the classical hyperbolic function, with half-maximal respiration rate observed at peak ADP concentrations close to 160-200 µM. At least two conclusions can be made from these results: 1) the relationship between cardiac tissue respiration rate and intracellular ADP concentration depends on how we calculate the latter (i.e., it is concept dependent) and, 2) clearly, ADP may participate in the feedback regulation of respiration in working cardiac cells, if its average or maximal peak values within the contraction cycle are taken into account. The second metabolic parameter that changes significantly is Pi (Fig. 9).


View larger version (15K):
[in this window]
[in a new window]
 
Fig. 6.   ADP profiles over cardiac cycle at different workloads. Profiles correspond to VO2 = 25 (low), 50 (medium), and 100 (high) µmol · g dry wt-1 · min-1. Profiles were obtained with 0-dimensional model.



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 7.   Relationship between maximal (max), minimal (min), and average (aver) levels of cytoplasmic ADP over a cardiac cycle and VO2 by working cardiac muscle. Considerable difference between maximal and average ADP levels is caused by significant change in ADP level during systole (see Fig. 11).



View larger version (14K):
[in this window]
[in a new window]
 
Fig. 8.   Relationship from Fig. 12 presented in more usual and conventional coordinates as relationship between VO2 by working cardiac muscle and maximal, minimal, and average levels of cytoplasmic ADP over a cardiac cycle. Considerable difference between maximal and average ADP levels is caused by significant change in ADP level during systole (see Fig. 11). Minimal ADP levels are close to those computed for CK equilibrium (see Fig. 6).



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 9.   Average computed ATP, Cr, PCr, and Pi levels over cardiac cycle as functions of VO2 by contracting heart muscle with increasing workload (see Fig. 7).

CK reaction is important for the function of normoxic heart at increased workloads. Two reactions of the intracellular energy transfer network (11, 32) are included in our model: CK and AK reactions. To evaluate their relative importance, we analyzed the situation when the activities of these enzymes were taken to be zero (modeling the experiments with inhibition of these enzymes or with their "knockout" due to genetic manipulations). With fully active CK, the blocking of AK activity did not alter the model behavior. This result confirms the conclusion of Dzeja et al. (11). Figure 4 shows the results of calculations with use of the model but the zero activity of CK and AK. This limited very much the workload and corresponding VO2 values that could be achieved by the heart: both stayed at only 20% of their maximally possible values. These values can be compared with those reported in three types of experiments: CK inhibition by iodoacetate (15), PCr replacement by feeding rats guanidinopropionate (39), or total knockout of CK (CK-deficient mice) (33). In all cases, the recorded maximal VO2 was 35-50 µmol · g dry wt-1 · min-1 and was significantly lower than the control values in normal hearts (if those were correctly studied at sufficiently high workloads). These experimental results are close to our calculated VO2 of 41 µmol · g dry wt-1 · min-1 (Fig. 4, open circle). Decreased work capacities of the skeletal muscles of CK-deficient mice have also been reported (36). Thus the model predicts that although the animals can easily survive without the MM isoform of CK or miCK activity, their work capacities are significantly compromised (15, 33, 36, 39).

Nonequilibrium state of CK and ADP concentration changes: importance of the intracellular diffusion rates. Figure 10 shows the numerical results of two- and zero-dimensional model calculations of the ADP concentration in intact cardiomyocytes, in myofibril core at 0.2 µm of the M line, where the activity of actomyosin ATPase is high (Fig. 1), and ADP levels at the boundary between the cytoplasm and the mitochondria during one cardiac contraction cycle. Calculations of these concentrations shown in Fig. 10 were made with a two-dimensional model for fixed energy fluxes equal to 100 µmol O2 · g dry wt-1 · min-1 by using the values of the diffusion coefficients characteristic of water. The dotted line in Fig. 10A was obtained for the case of infinitely fast diffusion of ADP (zero-dimensional model, see METHODS OF MODELING). In all cases, the ADP concentration is changed very significantly within the contraction cycle, its peak value exceeding about five times the resting (equilibrium) value of the ADP (see above). Remarkably, the two- and zero-dimensional models give close results, the maximal difference being only 10% of the peak value (solid and dotted lines in Fig. 10A). This means that significant variations of the ADP concentration in the cell are not related to problems of its diffusion but are intrinsic to the compartmentalized energy metabolism in the cell because of the rather low total CK activity and the specific intracellular distribution of CK. For the case of rapid diffusion of ADP (as in water), there is a small gradient of ADP between the myofibrillar core and the boundary of the myofibrillar space at the moment of its peak value (Fig. 10A).



View larger version (1313K):
[in this window]
[in a new window]
 
Fig. 10.   Time course of changes in ADP concentration in contracting cardiac cells. A: ADP concentrations [maximal in core of myofibrils (solid lines) and minimal at boundary between cytoplasm and mitochondria (dashed lines)] were calculated for 1 contraction cycle by use of normal diffusion coefficients (as in water), and average ADP concentrations (dotted lines) were obtained for infinitely fast diffusion in myofibrils [0-dimensional (0D) model]. B: same as A, but ADP levels were computed with 2-dimensional (2D) model with diffusion coefficients reduced 10-fold compared with those in water space.

Because there are some indications of a possible restriction of ADP diffusion in muscle cells (3), we repeated these calculations with a 10-fold decrease in the ADP diffusion coefficient (Fig. 10B). This increased the peak value of ADP concentration in the core of myofibrils and decreased it significantly at the boundary of the cytoplasm and the mitochondria; the ratio of ADP peak values was close to 3. Thus restriction of ADP diffusion increases its concentration gradients in the cells, as could be intuitively predicted.

The results in Fig. 10 were obtained for two points: for one point in the core of myofibrils, at 0.2 µm from the M line, and for another point at the mitochondrial surface. However, it was interesting to find out what may happen in other parts of the cell at the moment when the ADP concentration in the core of myofibrils is maximal. To answer this question, we fixed the time when the ADP concentration in the core of myofibrils reaches its maximal value, i.e., 40 ms, and analyzed the concentration gradients of ADP in the cell by using the two-dimensional model of compartmentalized energy transfer. These results are shown in Fig. 11. Figure 11A shows the ADP concentration profiles for normal values of its diffusion coefficient. The ADP concentrations at this moment in time are not equal in the cell and show concentration gradients. ADP concentration is maximal in the core of the sarcomere (0.2 µm at x-axis), where its value is 0.4 mM (Figs. 1 and 10) and decreases in the direction of the x- and y-axes, corresponding to directions to mitochondria and to the Z line. Solid and dashed lines in Fig. 10A correspond to maximal and minimal ADP concentrations in Fig. 11A. Figure 11B shows the contour plots of the ADP concentration in the myofibrillar space of cardiomyocytes for the case of restricted diffusion of ADP, when its diffusion coefficient was decreased 10-fold compared with its value in water. In accordance with the data from Fig. 10B, in this case, one observes very remarkable ADP concentration gradients within the cell.



View larger version (2735K):
[in this window]
[in a new window]
 
Fig. 11.   ADP distribution in 2 dimensions over myofibril and myoplasm in systole. Contour lines, concentrations (µM) at 40 ms after onset of contraction; arrows, direction and velocity of ADP flux. ADP concentrations are computed by 2-dimensional model with normal (water) diffusion coefficients (A) and with reduced (10 times less than normal) diffusion coefficients (B).

However, what happens with metabolites other than ADP that take part in the energy metabolism of the cells? The results of calculations, by use of the different reaction-diffusion models, of the ADP, ATP, Cr, PCr, and Pi concentrations in the core of myofibrils during the contraction cycle are shown in Fig. 12. The concentration changes within the contraction cycle are very significant only for ADP (~5 times), whereas at the same time, ATP concentration decreases only by 4%, Cr concentration increases 4%, and correspondingly PCr concentration decreases by 4% (Fig. 12). However, significant changes, by up to 25%, are seen for Pi concentration (Fig. 12) because of its low initial value. Assumption of an infinite rate of diffusion (zero-dimensional model) modifies only very slightly, within a range of 1%, these metabolite concentrations. Change of the dimension of the model has some influence only on the ADP concentration profile and has little influence on the others (Fig. 12). This means that the rate of intracellular diffusion of metabolites characteristic of the intracellular bulk water phase is sufficiently rapid for equilibration of the metabolites in the intracellular space, and further increase of diffusion coefficients does not alter significantly the metabolite profiles in the cells.


View larger version (25K):
[in this window]
[in a new window]
 
Fig. 12.   Dynamics of ADP, ATP, Cr, PCr, and Pi concentrations within cardiac cycle at myoplasmic side of mitochondrial outer membrane (at x = 0.5, y = 1.2 in Figs. 1 and 9) obtained from models with different spatial dimensions and diffusion coefficients. Solutions corresponding to 2-, 1-, and 0-dimensional models with normal (water) diffusion coefficients and 2- and 1-dimensional models with reduced (10 times less than normal) diffusion coefficients are presented.

Another question, however, relates to the state of the CK system. Significant variations in the myoplasmic ADP concentration within the contraction cycle shown in Figs. 10-13 demonstrate that the system is out of equilibrium. Figure 13 shows the two-dimensional model calculations of the metabolite profiles within a cardiac cycle for normal cellular activities of the compartmentalized CK in the heart and for the case when the CK activity was increased manyfold to ensure its equilibrium state. In the equilibrium, as expected, the levels of metabolites, including ADP, are practically constant, but in the real situation in the cell, the ADP concentration is close to its equilibrium only in the diastolic phase, its peak value exceeding the equilibrium five- to sixfold. There are no big differences between the real and equilibrium concentrations of Cr, PCr, and Pi; all these parameters are relatively stable, and the difference between real and equilibrium values is <10%. The concentration of ATP is constant in equilibrium and varies in the real situation <= 4.2%.


View larger version (23K):
[in this window]
[in a new window]
 
Fig. 13.   Dynamics of ADP, ATP, Cr, PCr, and Pi concentrations within cardiac cycle at myoplasmic side of mitochondrial outer membrane (at x = 0.5, y = 1.2 in Figs. 1 and 9) obtained by 2-dimensional models for increased (105 times compared with normal) MM CK activity (CK in equilibrium) and increased permeability of mitochondrial outer membrane (no outer membrane). For comparison, solutions of unmodified (normal) 2-dimensional model are also given. MM, MM isoform of CK.

There is, however, some influence of the permeability of the outer mitochondrial membrane for adenine nucleotides on the profiles of metabolites. The peak values of ADP are decreased by 20% after the outer mitochondrial membrane is "opened" (Fig. 13), and, correspondingly, the higher levels of PCr and lower levels of phosphate and free Cr are established (Fig. 13) within a cytoplasmic space.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

The results of this study show that the linear dependence of VO2 on workload of working cardiac muscle can be quantitatively described by the reaction-diffusion model of compartmentalized intracellular energy transfer, which is based on experimentally measured enzyme activities and accounts for intracellular localization of coupled CK. The model reproduces and explains the constant PCr-to-ATP ratios in a rather wide range of workloads. It shows that ATP, Cr, and PCr concentrations and their changes within the contraction cycle, calculated by using different modifications of the model, do not differ very much from those obtained on the basis of the classical assumption of CK equilibrium. However, the calculated myoplasmic ADP concentration shows very significant, workload-dependent variations within the cardiac cycle, and in the systolic phase its level may exceed its equilibrium value by an order of magnitude. In the metabolic research, the values of ATP, PCr, and Cr are the parameters determined experimentally (2, 12, 15, 17, 35, 38, 39). The results of this study show that interpretation of these data, i.e., calculation of ADP values, crucially depends on assumptions and models used. The validity of our new model is confirmed by quantitative description of experimental data for the whole physiological range of VO2 and, thus, energy fluxes in the heart cells in vivo. An analysis of the model shows that increasing infinitely the value of the diffusion constants does not alter significantly the character of changes of cytoplasmic concentrations of the metabolites within the cardiac cycle (Fig. 4). This means that in the bulk water phase the diffusion is already rapid and all calculated metabolite profiles in the given steady state are formed relative to the levels of the compartmentalized enzyme activities characteristic of the given cell. Because similar levels of ATP, Cr, and PCr are calculated by the equilibrium and compartmentalized energy transfer models, one may say that these metabolites have quasi-equilibrium steady-state levels (9). The basis of this phenomenon, according to our simulations, is that the diffusion of these compounds in the bulk water phase of the cytoplasm is sufficiently rapid and that, in the diastolic phase of the contraction cycle, the CK system in the myofibrillar space approaches an equilibrium. This allows one to use the system of ordinary differential equations (zero-dimensional model) in the place of the system of partial differential equations (two-dimensional model) to describe qualitatively these parameters of intracellular energy metabolism in the steady state. However, this conclusion of the small role of diffusion in the myoplasm in determining the metabolic profiles is based on the assumption that diffusion in the myoplasm may be modeled as a diffusion in the bulk water phase (26). Increasing experimental evidence shows that this may not be totally true in the living cells, where diffusion may be influenced by such factors as intracellular structures and the bound water structure (28; see below). However, the conclusion of the quasi-equilibrium is clearly not true for ADP concentration in the myoplasmic space. The level of ADP is lower than that of other metabolites by factor of 100, and its value changes manyfold in the contraction cycle with formation of the concentration gradients in the myofibrillar space revealed by two-dimensional modeling. Because the total content of adenine nucleotides is conserved in the model, ADP increases at the expense of ATP, but because of the 100-fold difference in initial levels, a decrease in ATP of only 3-4% is enough to change ADP by an order of magnitude. Our calculations of the diffusion profiles of ADP (Fig. 10) show that, because of its low concentration, ADP diffusion is not sufficiently rapid to equilibrate its concentration in the cells: there is a difference of 13% in ADP peak concentration between the core of myofibrils and the mitochondrial-myoplasmic border even for normal (characteristic for water) diffusion coefficients. This idea was first advocated by Kammermeier (20). Restriction of ADP diffusion results in formation of its very significant concentration gradients (Figs. 10 and 11). In this case, local ADP concentrations in myofibrils may reach very significant values (11).

Thus, in the nonresting state, where the energy fluxes are significant, the ADP concentration is clearly dissociated from quasi-equilibrium concentrations of other metabolites. It depends on the value of the energy flux in the system and on many other parameters (1).

Such a dissociation of ADP levels from levels of other metabolites that are close to the levels predicted by the CK equilibrium reopens the question of the importance of ADP in regulation of the rate of mitochondrial respiration. In the classical experiments of Balaban et al. (2), the increased VO2 values in hearts in vivo were observed at a constant PCr-to-ATP ratio. Very similar data were obtained in many other laboratories for heart and skeletal muscle (18, 35). The common conclusion from all these works was that there is no metabolic regulation of mitochondrial respiration by ADP in vivo, and a search for other alternative mechanisms, in particular by calcium ions, began. Our studies show that this conclusion, based on the assumption of CK equilibrium, is probably not valid because of the compartmentalization of the energy transfer system in the cell. In the compartmentalized CK system of the heart, the activities of different CK isozymes are very well fitted to generate the oscillations of ADP as a part of the metabolic feedback signal between contraction and energy production at the background of a constant level of other metabolites. These oscillations may be strongly amplified by the coupled reactions of aerobic PCr production in mitochondria (1). Furthermore, these oscillations of ADP concentration may be synchronized with localized changes of calcium in cytoplasm and in mitochondria, the latter modifying the activity of the enzymes of mitochondrial systems to increase the rate of ADP rephosphorylation (16), which in mitochondrial coupled reactions results in enhanced PCr production (31, 32). This explanation is in accord with the concept of parallel activation of energy-producing and -consuming processes proposed recently by Korzeniewski (21) also on the basis of the quantitative analysis. However, Korzeniewski and Froncisz (21, 22) completely omitted the CK reaction and, thus, the possible feedback metabolic regulation from their considerations. On the contrary, our calculations have been made for the constant maximal activity of mitochondrial oxidative phosphorylation. In fact, the phenomenon of parallel activation and any possible effects of calcium on respiration are ignored in our model. Without any need to account for these effects, our model explains quantitatively the workload dependencies of VO2 and remarkable metabolic stability, i.e., the unchanged levels of PCr for a significant range of work transitions (2, 38). However, the effects of calcium may be related to a change in the maximal velocity of respiration (thus activation of the respiration), but not to the problems of feedback regulation of mitochondrial oxidative phosphorylation in contracting cardiac cells. The real value of the rate of respiration is, according to our calculations, determined by the metabolic signal from cytoplasm. This is in accord with recent data of Brandes and Bers (4). This signal, in the form of changing ADP and Pi concentrations, amplified in the coupled CK reaction (32), determines the fraction of maximal reaction velocity used in any metabolic conditions. Thus, taken together, these two theories may explain the feedback regulation of the respiration under conditions of metabolic stability by synchronized oscillations of the calcium and ADP concentrations. These concentration changes may well be localized in the intracellular microcompartments. As shown recently by Rizzuto et al. (29), the calcium concentration changes may be very significantly localized in the mitochondrial and surrounding microcompartments because of close contacts between mitochondria and sarcoplasmic reticulum. Our results and conclusions are also in accord with the recent discovery of subcellular metabolic transients and mitochondrial redox waves in cardiomyocytes by O'Rourke et al. (27) and Romashko et al. (30). Using confocal imaging of flavoprotein redox potential and mitochondrial membrane potential, they showed that substrate deprivation leads to subcellular heterogeneity of mitochondrial energization in the cell and propagation of a redox wave at ~2 µm/s, which is comparable to the average diffusion rate of small molecules such as ADP, if the diffusion coefficient is decreased by a factor of 5 compared with that for water (30). The authors concluded that intracellular control of mitochondrial function involves diffusible cytoplasmic messengers. The role of calcium as a mediator was excluded, since the redox oscillations were independent of calcium concentration (30). However, rapid release of ADP by flash photolysis of intracellular caged ADP induced metabolic oscillations (27). It has been proposed that cytoplasmic structures may control mitochondria and that there exists an organized network of energy transfer and feedback signal transduction in the cells (11, 32) by a mechanism of "vectorial ligand conduction through organized cytoplasmic multienzyme systems" [a definition proposed by Mitchell (24)], including CK and AK and probably glycolytic and other systems. In this system, free diffusion of ADP is replaced by its vectorial conduction, along with other ligands, by enzymes, probably organized in association with the cytoskeleton. In fact, free diffusion of ADP itself seems to be very limited in the muscle cells. Multiple experiments carried out by Ventura-Clapier and by Bessman et al. (see Ref. 32 for references) on the skinned muscle fibers showed that relaxation from the rigor state was more effective by an order of magnitude by PCr because of local rephosphorylation of ADP than by externally added ATP. This shows that, in the absence of PCr, ADP is accumulating in the vicinity of the myosin-active centers and slows the cross-bridge detachment, and diffusion of ATP into fibers is not as effective as diffusion of PCr. Bereiter-Hahn and Voth (3) described a very interesting phenomenon of delayed morphological response of mitochondria to ATP or ADP microinjected into living cells, indicating that diffusion of these molecules is ~10 times slower than that of other small hydrophilic molecules. This kind of structural organization of the cytoplasmic space in the cells in vivo is not yet included in our model. Obviously, this kind of organization localizes the changes of the ADP concentration in the cytoplasmic space (ADP "sparks"), which may propagate as a metabolic wave, giving rise to oscillations of other metabolites (ligands) (8, 29, 32). Clearly, the next very interesting step is to include this kind of organized and microcompartmentalized vectorial process in our model of energy transfer and feedback regulation of respiration. This may help in analysis of the recently observed metabolic transients and redox waves in cardiac cells (27, 30).

In conclusion, the new version of the mathematical model of compartmentalized energy transfer in cardiac cells describes quantitatively the in vitro data on isolated mitochondria and the in vivo metabolic data on working perfused rat heart. It describes 1) dependence of the rate of PCr production in the coupled CK reaction on oxidative phosphorylation in vitro, 2) linear dependence of VO2 by perfused heart on the workload, 3) reduction of the maximal workload of perfused heart after inhibition of CK, and 4) experimentally observed metabolic stability (constant PCr-to-ATP ratio) of heart cells at different workloads. That means that a feedback mechanism can regulate the rate of oxidative phosphorylation in a wide range of workloads and keep the levels of metabolites, i.e., Cr, PCr, and ATP, at metastable steady-state levels and, thus, can explain the apparent controversy widely discussed in the literature (2, 4, 15-23, 30, 32-35, 38). The model shows that the amplitude of ADP oscillations due to the nonequilibrium state of the CK system increases in cytoplasm with a decrease in the diffusion coefficient and that ADP gradients are larger along than across the myofibrils. The latter conclusions point to the potential importance of developing independent experimental methods for determination of the ADP concentrations in different cellular compartments, which is still a challenge in cellular bioenergetics.

However, to evaluate quantitatively the contribution of changes in cytoplasmic ADP and Pi concentrations to regulation of the rate of mitochondrial respiration in vivo, further analysis of the model is required by using approaches similar to those used in metabolic control analysis.


    APPENDIX
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

Here we describe the mathematical models. Depending on dimensions, we call these two- and zero-dimensional models. The two-dimensional model is a continuous version of the discrete one-dimensional model (1) obtained by adding the longitudinal dimension. This is a system of partial differential equations involving two spatial coordinates and time. With the assumption of spatial homogeneity of the metabolites within the myocyte, i.e., infinitely fast diffusion, the two-dimensional model is simplified to a system of ordinary differential equations involving only time dependence, hence, the name zero-dimensional model. The values of the parameters used in the models are listed in Table 1.

Two-dimensional model. The dynamics of the metabolites in the myofibril and myoplasm are described by the following system of reaction-diffusion equations
<FR><NU>∂ATP</NU><DE>∂<IT>t</IT></DE></FR> = <IT>D</IT><SUB>ATP</SUB>∇<SUP>2</SUP>ATP − <IT>G</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> + <IT>G</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> − <IT>G</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) (1)

<FR><NU>∂ADP</NU><DE>∂<IT>t</IT></DE></FR> = <IT>D</IT><SUB>ADP</SUB>∇<SUP>2</SUP>ADP + <IT>G</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> − 2<IT>G</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> + <IT>G</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) (2)

<FR><NU>∂AMP</NU><DE>∂<IT>t</IT></DE></FR> = <IT> 
D</IT><SUB>AMP</SUB>∇<SUP>2</SUP>AMP + <IT>G</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> (3)

<FR><NU>∂PCr</NU><DE>∂<IT>t</IT></DE></FR> = <IT>D</IT><SUB>PCr</SUB>∇<SUP>2</SUP>PCr + <IT>G</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> (4)

<FR><NU>∂Cr</NU><DE>∂<IT>t</IT></DE></FR> = <IT>D</IT><SUB>Cr</SUB>∇<SUP>2</SUP>Cr − <IT>G</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> (5)

<FR><NU>∂P<SUB>i</SUB></NU><DE>∂<IT>t</IT></DE></FR> = <IT>D</IT><SUB>P<SUB>i</SUB></SUB>∇<SUP>2</SUP>P<SUB>i</SUB> + <IT>G</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) (6)
where ATP, ADP, AMP, PCr, Cr, and Pi are the total concentrations of the metabolites in the myofibril and myoplasm, DMet is the diffusion coefficient of the metabolite (Met, i.e., ATP, ADP, AMP, PCr, Cr, Pi), vCK is the CK reaction rate, vAK is the AK reaction rate, and HATP is the rate of ATP hydrolysis in the myofibril. The factors GCK, GAK, and GH describe the spatial inhomogeneities of myofibrillar CK, AK, and ATPase, respectively. The symbol nabla 2 denotes the Laplace operator.

The vCK is described by Eq. 7
<IT>v</IT><SUB>CK</SUB> = <FENCE><IT>V</IT><SUB>1</SUB> <FR><NU>mATP ⋅ Cr</NU><DE><IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB></DE></FR> − <IT>V</IT><SUB>−1</SUB> <FR><NU>mADP ⋅ PCr</NU><DE><IT>K</IT><SUB>ic</SUB><IT>K</IT><SUB>d</SUB></DE></FR></FENCE> &cjs0823;   DenCK (7)
where
DenCK = 1 + <FR><NU>Cr</NU><DE><IT>K</IT><SUB>ib</SUB></DE></FR> + <FR><NU>PCr</NU><DE><IT>K</IT><SUB>id</SUB></DE></FR> + mATP <FENCE><FR><NU>1</NU><DE><IT>K</IT><SUB>ia</SUB></DE></FR> + <FR><NU>Cr</NU><DE><IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB></DE></FR></FENCE>

+ mADP <FENCE><FR><NU>1</NU><DE><IT>K</IT><SUB>ic</SUB></DE></FR> + <FR><NU>PCr</NU><DE><IT>K</IT><SUB>id</SUB><IT>K</IT><SUB>c</SUB></DE></FR> + <FR><NU>Cr</NU><DE><IT>K</IT><SUB>ic</SUB><IT>K</IT><SUB>Ib</SUB></DE></FR></FENCE> (8)
and
<IT>K</IT><SUB>c</SUB> = <IT>K</IT><SUB>ic</SUB><IT>K</IT><SUB>d</SUB>/<IT>K</IT><SUB>id</SUB> (9)
where mATP and mADP are the concentrations of the Mg2+-bound forms of ATP and ADP, respectively, and fATP and fADP are the Mg2+-free forms. The total concentrations of adenylates and their Mg2+-bound forms are related as follows
fATP = <FR><NU>ATP</NU><DE>1 + Mg/<IT>K</IT><SUB>DT</SUB></DE></FR> (10)

mATP = ATP − fATP (11)

fADP = <FR><NU>ADP</NU><DE>1 + Mg/<IT>K</IT><SUB>DD</SUB></DE></FR> (12)

mADP = ADP − fADP (13)
where Mg is the level of free Mg2+ and KDT and KDD are Mg2+ dissociation constants. The level of free Mg2+ and its dissociation constant for the matrix differ from their corresponding values for the rest of the cell (22). The vAK is described as follows
<IT>v</IT><SUB>AK</SUB> = <IT>k</IT><SUB>fa</SUB> ⋅ fADP ⋅ mADP − <IT>k</IT><SUB>ba</SUB> ⋅ mATP ⋅ AMP (14)
where kfa and kba are the forward and backward reaction rate constants for the AK reaction.

We approximate the hydrolysis rate by a piecewise linear time-periodic function HATP. The rate increases from 0 to HATP(max) during the first 30 ms, decreases to 0 during the next 30 ms, and remains 0 until the end of the cardiac cycle at 180 ms
<IT>H</IT><SUB>ATP</SUB> = <FENCE> <AR><R><C><FR><NU>t</NU><DE>30.0</DE></FR> <IT>H</IT><SUB>ATP(max)</SUB>,    if 0 ≤ <IT>t</IT> < 30</C></R><R><C><FR><NU>60.0 − <IT>t</IT></NU><DE>30.0</DE></FR> <IT>H</IT><SUB>ATP(max)</SUB>,  if 30 ≤ <IT>t</IT> < 60</C></R><R><C>0.0,        if 60 ≤ <IT>t</IT> < 180</C></R></AR></FENCE> (15)
where t is elapsed time from the beginning of the cycle (in ms).

The reaction rates of CK, AK, and ATPase are proportional to the concentrations of the corresponding enzymes. The spatial distribution of CK (37) is described by
<IT>G</IT><SUB>CK</SUB>(<IT>x</IT>, <IT>y</IT>) 

= <FENCE><AR><R><C>0.94<FENCE>1 − <FR><NU><IT>x</IT></NU><DE>0.375</DE></FR></FENCE>, if <IT>x</IT> < 0.375 and <IT>y</IT> < 1</C></R><R><C>1.61 <FR><NU><IT>x</IT> − 0.375</NU><DE>0.375</DE></FR>, if 0.375 ≤ <IT>x</IT> < 0.75 and <IT>y</IT> < 1</C></R><R><C>1.61,      if <IT>x</IT> ≥ 0.75 or <IT>y</IT> ≥ 1</C></R></AR></FENCE> (16)
the distribution of AK by
<IT>G</IT><SUB>AK</SUB>(<IT>x</IT>, <IT>y</IT>) = <IT>G</IT><SUB>CK</SUB>(<IT>x</IT>, <IT>y</IT>) (17)
and the distribution of ATPase (5) by
<IT>G</IT><SUB>H</SUB>(<IT>x</IT>, <IT>y</IT>)

= <FENCE> <AR><R><C><FR><NU>4</NU><DE>3</DE></FR> <FR><NU><IT>x</IT></NU><DE>0.125</DE></FR> , if <IT>x</IT> < 0.125 and <IT>y</IT> < 1</C></R><R><C><FR><NU>4</NU><DE>3</DE></FR> ,if 0.125 ≤ <IT>x</IT> < 0.75 and <IT>y</IT> < 1</C></R><R><C><FR><NU>4</NU><DE>3</DE></FR> <FENCE>1 − <FR><NU><IT>x</IT> − 0.75</NU><DE>0.125</DE></FR></FENCE> , if 0.75 ≤ <IT>x</IT> < 0.875 and <IT>y</IT> < 1</C></R><R><C>0,if <IT>x</IT> ≥ 0.875 or <IT>y</IT> ≥ 1</C></R></AR></FENCE> (18)
where x and y are spatial coordinates (in µm): x corresponds to the longitudinal direction (increases from the M line to the Z line), and y corresponds to the transverse direction (increases from the core of myofibrils to the mitochondrion); the point with coordinates (0,0) is located at the core of myofibrils and the M line. The mean values of GCK and GAK over the myofibril and myoplasm (0 <=  x <=  1, 0 <=  y <=  1.2) and GH over the myoplasm (0 <=  x <=  1, 0 <=  y <=  1) are equal to unity.

Because of the symmetry of the problem, we use no flux condition at the boundaries x = 0, x = 1, and y = 0. The restricted diffusion through the mitochondrial outer membrane takes place at the boundary y = 1.2 and is calculated by
<IT>D</IT><SUB>ATP</SUB> <FR><NU>∂ATP</NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>ATP</SUB>(ATP<SUB>i</SUB> − ATP) (19)

<IT>D</IT><SUB>ADP</SUB> <FR><NU>∂ADP</NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>ADP</SUB>(ADP<SUB>i</SUB> − ADP) (20)

<IT>D</IT><SUB>AMP</SUB> <FR><NU>∂AMP</NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>AMP</SUB>(AMP<SUB>i</SUB> − AMP) (21)

<IT>D</IT><SUB>PCr</SUB> <FR><NU>∂PCr</NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>PCr</SUB>(PCr<SUB>i</SUB> − PCr) (22)

<IT>D</IT><SUB>Cr</SUB> <FR><NU>∂Cr</NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>Cr</SUB>(Cr<SUB>i</SUB> − Cr) (23)

<IT>D</IT><SUB>P<SUB>i</SUB></SUB> <FR><NU>∂P<SUB>i</SUB></NU><DE>∂<B><A><AC>n</AC><AC>&cjs1164;</AC></A></B></DE></FR> = <IT>R</IT><SUB>P<SUB>i</SUB></SUB>(P<SUB>i,i</SUB> − P<SUB>i</SUB>) (24)
where Meti is the concentration of the metabolite in the mitochondrial IM space and RMet is the permeability of the mitochondrial outer membrane.

The concentrations of the metabolites in the IM space are changed because of the diffusion through the outer membrane, CK and AK reactions, ATP/ADP translocation, and Pi transportation through the inner membrane according to the following equations
<FR><NU>dATP<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>ATP</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (ATP − ATP<SUB>i</SUB>) + <IT>v</IT><SUB>ANT</SUB>/FV<SUB>i</SUB>

 − <IT>v</IT><SUB>MiCK</SUB>/FV<SUB>i</SUB> + <IT>v</IT><SUB>AKmit</SUB>/FV<SUB>i</SUB> (25)

<FR><NU>dADP<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>ADP</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (ADP − ADP<SUB>i</SUB>) − 2<IT>v</IT><SUB>AKmit</SUB>/FV<SUB>i</SUB>

 − <IT>v</IT><SUB>ANT</SUB>/FV<SUB>i</SUB> + <IT>v</IT><SUB>MiCK</SUB>/FV<SUB>i</SUB> (26)

<FR><NU>dAMP<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>AMP</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (AMP − AMP<SUB>i</SUB>) + <IT>v</IT><SUB>AKmit</SUB>/FV<SUB>i</SUB> (27)

<FR><NU>dPCr<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>PCr</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (PCr − PCR<SUB>i</SUB>) + <IT>v</IT><SUB>AKmit</SUB>/FV<SUB>i</SUB> (28)

<FR><NU>dCr<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>Cr</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (Cr − Cr<SUB>i</SUB>) − <IT>v</IT><SUB>MiCK</SUB>/FV<SUB>i</SUB> (29)

<FR><NU>dP<SUB>i,i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>R</IT><SUB>Pi</SUB></NU><DE><IT>L</IT><SUB>i</SUB></DE></FR> (P<SUB>i</SUB> − P<SUB>i,i</SUB>) − <IT>v</IT><SUB>PI</SUB>/FV<SUB>i</SUB> (30)
where vMiCK is MiCK reaction rate, vANT is net rate of ATP export by ANT from the matrix, vAKmit is AK reaction rate (see Eq. 14), vPi is rate of the phosphate carrier in the inner membrane, Li is width of the layer between the inner and outer membranes, and FVi is fractional volume of the IM space with respect to the cell volume.

The MiCK reaction rate (vMiCK) is divided into vMiCK,G, which involves microcompartment ATP(ADP), and vMiCK,I, which involves ATP(ADP), from the IM space
<IT>v</IT><SUB>MiCK</SUB> = <IT>v</IT><SUB>MiCK,G</SUB> + <IT>v</IT><SUB>MiCK,I</SUB> (31)
The reaction rates are as follows
<IT>v</IT><SUB>MiCK,G</SUB> = <FENCE><IT>V</IT><SUB>1</SUB> <FR><NU>mATP<SUB>g</SUB> ⋅ Cr<SUB>i</SUB></NU><DE><IT>K</IT><SUB>ia,G</SUB><IT>K</IT><SUB>b</SUB></DE></FR> − <IT>V</IT><SUB>−1</SUB> <FR><NU>mADP<SUB>g</SUB> ⋅ PCr<SUB>i</SUB></NU><DE><IT>K</IT><SUB>ic,G</SUB><IT>K</IT><SUB>d</SUB></DE></FR></FENCE> &cjs0823;   DenMiCK (32)

<IT>v</IT><SUB>MiCK,I</SUB> = <FENCE><IT>V</IT><SUB>1</SUB> <FR><NU>mATP<SUB>i</SUB> ⋅ Cr<SUB>i</SUB></NU><DE><IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB></DE></FR> − <IT>V</IT><SUB>−1</SUB> <FR><NU>mADP<SUB>i</SUB> ⋅ PCr<SUB>i</SUB></NU><DE><IT>K</IT><SUB>ic</SUB><IT>K</IT><SUB>d</SUB></DE></FR></FENCE> &cjs0823;   DenMiCK (33)
where


DenMiCK = 1 + Cr<SUB>i</SUB>/<IT>K</IT><SUB>ib</SUB> + PCr<SUB>i</SUB>/<IT>K</IT><SUB>id</SUB> + (mATP<SUB>g</SUB>/<IT>K</IT><SUB>ia,G</SUB> + mATP<SUB>i</SUB>/<IT>K</IT><SUB>ia</SUB>) ⋅ (1 + Cr<SUB>i</SUB>/<IT>K</IT><SUB>b</SUB>) 

+ (mADP<SUB>g</SUB>/<IT>K</IT><SUB>ic,G</SUB> + mADP<SUB>i</SUB>/<IT>K</IT><SUB>ic</SUB>) ⋅ (1 + PCr<SUB>i</SUB>/<IT>K</IT><SUB>d</SUB> + Cr<SUB>i</SUB>/<IT>K</IT><SUB>Ib</SUB>) (34)

Restricted diffusion between the microcompartment and the IM space is governed by
<IT>v</IT><SUB>exchATP</SUB> = <IT>R</IT><SUB>exchATP</SUB> ⋅ (ATP<SUB>i</SUB> − ATP<SUB>g</SUB>) (35)
for ATP and by
<IT>v</IT><SUB>exchADP</SUB> = <IT>R</IT><SUB>exchADP</SUB> ⋅ (ADP<SUB>i</SUB> − ADP<SUB>g</SUB>) (36)
for ADP.

For the oxidative phosphorylation, we use the model by Korzeniewski (21)
<FR><NU>dUQ</NU><DE>d<IT>t</IT></DE></FR> = (<IT>v</IT><SUB>C3</SUB> − <IT>v</IT><SUB>C1</SUB>)/FV<SUB>x</SUB> (37)

<FR><NU>dc<SUP>3+</SUP></NU><DE>d<IT>t</IT></DE></FR> = 2(2<IT>v</IT><SUB>O<SUB>2</SUB></SUB> − <IT>v</IT><SUB>C3</SUB>)/FV<SUB>x</SUB> (38)

<FR><NU>dNAD<SUP>+</SUP></NU><DE>d<IT>t</IT></DE></FR> = (<IT>v</IT><SUB>C1</SUB> − <IT>v</IT><SUB>dh</SUB>)/FV<SUB>x</SUB> (39)

<FR><NU>dH<SUB>x</SUB></NU><DE>d<IT>t</IT></DE></FR> = −[2(2 + 2<IT>u</IT>)<IT>v</IT><SUB>O<SUB>2</SUB></SUB> + 4<IT>v</IT><SUB>C1</SUB> + (4 − 2<IT>u</IT>)<IT>v</IT><SUB>C3</SUB> − <IT>n</IT><SUB>a</SUB><IT>v</IT><SUB>sn</SUB> − +<IT>v</IT><SUB>ANT</SUB><IT>u</IT> − (1 − <IT>u</IT>)<IT>v</IT><SUB>PI</SUB> − <IT>v</IT><SUB>leak</SUB>]/(FV<SUB>x</SUB> ⋅ <IT>r</IT><SUB>buff</SUB>) (40)

<FR><NU>dATP<SUB>x</SUB></NU><DE>d<IT>t</IT></DE></FR> = (<IT>v</IT><SUB>sn</SUB> − <IT>v</IT><SUB>ANT</SUB>)/FV<SUB>x</SUB> (41)

<FR><NU>dP<SUB>i,x</SUB></NU><DE>d<IT>t</IT></DE></FR> = (<IT>v</IT><SUB>PI</SUB> − <IT>v</IT><SUB>sn</SUB>)/FV<SUB>x</SUB> (42)
where UQ is the oxidized form of coenzyme Q, c3+ is the oxidized form of cytochrome c, NAD+ is NAD+ in the matrix, Hx, ATPx, and Pi,x are H+, ATP, and Pi in the matrix, and FVx is fractional volume of the mitochondrial matrix with respect to the cell volume. The following reaction rates are used in the above equations (21)
<IT>v</IT><SUB>dh</SUB> = <FR><NU><IT>k</IT><SUB>dh</SUB></NU><DE><FENCE>1 + <FR><NU><IT>k</IT><SUB>mN</SUB></NU><DE>NAD<SUP>+</SUP>/NADH</DE></FR></FENCE><SUP> <IT>p</IT><SUB>D</SUB></SUP></DE></FR> (43)
for substrate dehydrogenation (dh)
<IT>v</IT><SUB>C1</SUB> = <IT>k</IT><SUB>C1</SUB>(<IT>E</IT><SUB>u</SUB> − <IT>E</IT><SUB>n</SUB> − 2&Dgr;<IT>p</IT>) (44)
for complex I (C1)
<IT>v</IT><SUB>C3</SUB> = <IT>k</IT><SUB>C3</SUB>[<IT>E<SUB>c</SUB> − E</IT><SUB>u</SUB> − (2 − <IT>u</IT>)&Dgr; <IT>p</IT>] (45)
for complex III (C3)
<IT>v</IT><SUB>O<SUB>2</SUB></SUB> = <FR><NU>1</NU><DE>2</DE></FR> <IT>k</IT><SUB>C4</SUB> ⋅ a<SUP>2+</SUP> ⋅ c<SUP>2+</SUP> <FR><NU>O<SUB>2</SUB></NU><DE>O<SUB>2</SUB> + <IT>k</IT><SUB>mO</SUB></DE></FR> (46)
for complex IV (O2)
<IT>v</IT><SUB>sn</SUB> = <IT>k</IT><SUB>sn</SUB> <FR><NU>10<SUP>&Dgr;<IT>G</IT><SUB>sn</SUB>/<IT>Z</IT></SUP> − 1</NU><DE>10<SUP>&Dgr;<IT>G</IT><SUB>sn</SUB>/<IT>Z</IT></SUP> + 1</DE></FR> (47)
for ATP synthase (sn)
<IT>v</IT><SUB>PI</SUB> = <IT>v</IT><SUB>f,PI</SUB>H<SUB>ext</SUB> <FR><NU><IT>P</IT><SUB>i,i</SUB></NU><DE>1 + 10<SUP>pH<SUB>ext</SUB> − p<IT>K</IT><SUB>a</SUB></SUP></DE></FR> − <IT>v</IT><SUB>b,PI</SUB>H<SUB>x</SUB> <FR><NU>P<SUB>i,x</SUB></NU><DE>1 + 10<SUP>pH<SUB>x</SUB> − p<IT>K</IT><SUB>a</SUB></SUP></DE></FR> (48)
for phosphate carrier (PI), and
<IT>v</IT><SUB>leak</SUB> = <IT>k</IT><SUB>L1</SUB>(10<SUP><IT>k</IT><SUB>L2</SUB>&Dgr;<IT>p</IT></SUP> − 1) (49)
for proton leak (leak). In Eqs. 43-49, we used the total (tot) metabolite pools for substrate (NAD), coenzyme Q (UQ), cytochrome c (c), cytochrome a3 (a), and adenine nucleotides (in matrix, Ax) as follows
NAD<SUB>tot</SUB> = NAD<SUP>+</SUP> + NADH (50)

UQ<SUB>tot</SUB> = UQ + UQH<SUB>2</SUB> (51)

c<SUB>tot</SUB> = c<SUP>3+</SUP> + c<SUP>2+</SUP> (52)

a<SUB>tot</SUB> = a<SUP>3/2</SUP> + a<SUP>2+</SUP> (53)

A<SUB>x</SUB> = ATP<SUB>x</SUB> + ADP<SUB>x</SUB> (54)
The concentrations of Mg2+-bound and Mg2+-free forms of adenine nucleotides in the matrix are obtained using Eqs. 10-13. For computing pH in the matrix (pHx), pH difference across the inner membrane (Delta pH), protonmotive force (Delta p), electrical potential (Delta psi ), and buffering rate (rbuff), the following equations were used
pH<SUB>x</SUB> = −log(&eegr; ⋅ H<SUB>x</SUB>) (55)

&Dgr;pH = <IT>Z</IT>(pH<SUB>ext</SUB> − pH<SUB>x</SUB>) (56)

&Dgr;<IT>p</IT> = <FR><NU>1</NU><DE>1 − <IT>u</IT></DE></FR> &Dgr;pH (57)

&Dgr;&PSgr; = −(&Dgr;<IT>p</IT> − &Dgr;pH) (58)

H<SUB>ext</SUB> = &eegr; ⋅ 10<SUP>−pH<SUB>ext</SUB></SUP> (59)

<IT>r</IT><SUB>buff</SUB> = <IT>r</IT><SUB>buff,0</SUB> ⋅ &eegr; ⋅ <FR><NU>dpH<SUB>x</SUB></NU><DE>dH<SUB>x</SUB></DE></FR> (60)
where the coefficient u = Delta Psi /Delta p is kept constant. Because of the high permeability of the mitochondrial outer membrane for protons, pH in IM space (pHext) is considered to be equal to that in the cytoplasm and is considered to be constant. The coefficient Z is used to convert between molar concentrations and the concentrations used in the model (µM). The coeffiient Z is expressed as
<IT>Z</IT> = ln(10) <FR><NU><IT>RT</IT></NU><DE><IT>F</IT></DE></FR> (61)
where R is the gas constant, T is the absolute temperature, and F is Faraday's number. Throughout this study, the notations ln and log represent logarithms on bases e = 2.718,... and 10, respectively.

Description of ANT kinetics is based on the phenomenological model by Korzeniewski (21). The net rate of ATP export by ANT from the matrix to the microcompartment and the IM space is
<IT>v</IT><SUB>ANT</SUB> = <IT>v</IT><SUP>X→G</SUP><SUB>ANT</SUB> + <IT>v</IT><SUP>X→I</SUP><SUB>ANT</SUB> (62)
where vright-arrow  GANT is the rate of ATP export by ANT from the matrix to the microcompartment and vright-arrow  IANT is the rate of ATP export by ANT from the matrix directly to the IM space. These rates are described by the following kinetic equations
<IT>v</IT><SUP>X→G</SUP><SUB>ANT</SUB> = <IT>V</IT><SUB>ANT</SUB> ⋅ <FR><NU>fADP<SUB>g</SUB></NU><DE><IT>K</IT><SUB>g</SUB>(1 + fADP<SUB>g</SUB>/<IT>K</IT><SUB>g</SUB> + fADP<SUB>i</SUB>/<IT>K</IT><SUB>i</SUB>)</DE></FR> × <FENCE><FR><NU>fADP<SUB>g</SUB></NU><DE>fATP<SUB>g</SUB>10<SUP>−0.35&Dgr;&PSgr;/<IT>Z</IT></SUP> + fADP<SUB>g</SUB></DE></FR> − <FR><NU>fADP<SUB>x</SUB></NU><DE>fATP<SUB>x</SUB>10<SUP>−0.65&Dgr;&PSgr;/<IT>Z</IT></SUP> + fADP<SUB>x</SUB></DE></FR></FENCE> (63)

<IT>v</IT><SUP>X→I</SUP><SUB>ANT</SUB> = <IT>V</IT><SUB>ANT</SUB> ⋅ <FR><NU>fADP<SUB>i</SUB></NU><DE><IT>K</IT><SUB>i</SUB>(1 + fADP<SUB>g</SUB>/<IT>K</IT><SUB>g</SUB> + fADP<SUB>i</SUB>/<IT>K</IT><SUB>i</SUB>)</DE></FR> × <FENCE><FR><NU>fADP<SUB>i</SUB></NU><DE>fATP<SUB>i</SUB>10<SUP>−0.35&Dgr;&PSgr;/<IT>Z</IT></SUP> + fADP<SUB>i</SUB></DE></FR> − <FR><NU>fADP<SUB>x</SUB></NU><DE>fATP<SUB>x</SUB>10<SUP>−0.65&Dgr;&PSgr;/<IT>Z</IT></SUP> + fADP<SUB>x</SUB></DE></FR></FENCE> (64)

The capacity of the microcompartment is considered to be infinitesimal. Therefore, for ATP and ADP, the influx of the metabolite is balanced by its efflux
−<IT>v</IT><SUB>MiCK,G</SUB> + <IT>v</IT><SUP>X→G</SUP><SUB>ANT</SUB> + <IT>v</IT><SUB>exchATP</SUB> = 0 (65)

<IT>v</IT><SUB>MiCK,G</SUB> + <IT>v</IT><SUP>X→G</SUP><SUB>ANT</SUB> + <IT>v</IT><SUB>exchADP</SUB> = 0 (66)
Equations 65 and 66 can be solved for ATPg and ADPg by using, for example, the Newton-Raphson method.

The thermodynamic span of ATP synthase is calculated as follows
&Dgr;<IT>G</IT><SUB>sn</SUB> = <IT>n</IT><SUB>a</SUB>&Dgr;<IT>p</IT> − <FR><NU>&Dgr;<IT>G</IT><SUB>0</SUB></NU><DE><IT>F</IT></DE></FR> + <IT>Z</IT> log <FR><NU>&eegr;ATP<SUB>x</SUB></NU><DE>ADP<SUB>x</SUB>P<SUB>i,x</SUB></DE></FR> (67)
where na and Delta G0 are the H+/ATP stoichiometry and standard free energy for ATP synthase, respectively.

In Eqs. 44 and 45, the NAD, coenzyme Q, and cytochrome c redox potentials are used
<IT>E</IT><SUB>n</SUB> = <IT>E</IT><SUB>n,0</SUB> + <FR><NU><IT>Z</IT></NU><DE>2</DE></FR> log <FR><NU>NAD<SUP>+</SUP></NU><DE>NADH</DE></FR> (68)

<IT>E</IT><SUB>u</SUB> = <IT>E</IT><SUB>u,0</SUB> + <FR><NU><IT>Z</IT></NU><DE>2</DE></FR> log <FR><NU>UQ</NU><DE>UQH<SUB>2</SUB></DE></FR> (69)

<IT>E</IT><SUB>c</SUB> = <IT>E</IT><SUB>c,0</SUB> + <FR><NU><IT>Z</IT></NU><DE>2</DE></FR> log <FR><NU>c<SUP>3+</SUP></NU><DE>c<SUP>2+</SUP></DE></FR> (70)
respectively. For computing the concentration of reduced cytochrome a3 used in Eq. 46, the redox potential of cytochrome a3
<IT>E</IT><SUB>a</SUB> = <IT>E</IT><SUB>c</SUB> + (1 + <IT>u</IT>)&Dgr;<IT>p</IT> (71)
is substituted into the relation
<IT>a</IT><SUP>3/2</SUP> = 10<SUP>(<IT>E</IT><SUB>a</SUB> − <IT>E</IT><SUB>a,0</SUB>)<IT>Z</IT></SUP> (72)
which gives the result (together with Eq. 53).

With the assumption of infinitely fast diffusion, the two-dimensional model can be simplified to the zero-dimensional model. In this model the concentrations of the metabolites are spatially homogeneous and all the partial differential equations can be reduced to ordinary differential equations, hence, the name. The zero-dimensional model is obtained by setting
∇<SUP>2</SUP>Met = 0 (73)
in Eqs. 1-6, integrating them over the region of myofibril and myoplasm, and dividing by the area of this region. The integration of space-dependent functions GCK, GAK, and GH leads us to
<IT>F</IT><SUB>CK</SUB> = <FR><NU>∫∫ <IT>G</IT><SUB>CK</SUB>(<IT>x</IT>, <IT>y</IT>) d<IT>x</IT> d<IT>y</IT></NU><DE>∫∫ d<IT>x</IT> d<IT>y</IT></DE></FR> = 1 (74)

<IT>F</IT><SUB>AK</SUB> = <FR><NU>∫∫ <IT>G</IT><SUB>AK</SUB>(<IT>x</IT>, <IT>y</IT>) d<IT>x</IT> d<IT>y</IT></NU><DE>∫∫ d<IT>x</IT> d<IT>y</IT></DE></FR> = 1 (75)
and
<IT>F</IT><SUB>H</SUB> = <FR><NU>∫∫ <IT>G</IT><SUB>H</SUB>(<IT>x</IT>, <IT>y</IT>) d<IT>x</IT> d<IT>y</IT></NU><DE>∫∫ d<IT>x</IT> d<IT>y</IT></DE></FR> = <FR><NU>1</NU><DE>1.2</DE></FR> (76)

Zero-dimensional model. The concentrations of the free metabolites in cardiac myocytes are described by the following equations
<FR><NU>dATP</NU><DE>d<IT>t</IT></DE></FR> = −<IT>F</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> + <IT>F</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> (77)

 − <IT>F</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) + <FR><NU><IT>R</IT><SUB>ATP</SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (ATP<SUB>i</SUB> − ATP)

<FR><NU>dADP</NU><DE>d<IT>t</IT></DE></FR> = <IT>F</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> − 2<IT>F</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> (78)

 + <IT>F</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) + <FR><NU><IT>R</IT><SUB>ADP</SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (ADP<SUB>i</SUB> − ADP)

<FR><NU>dAMP</NU><DE>d<IT>t</IT></DE></FR> = <IT>F</IT><SUB>AK</SUB><IT>v</IT><SUB>AK</SUB> + <FR><NU><IT>R</IT><SUB>AMP</SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (AMP<SUB>i</SUB> − AMP) (79)

<FR><NU>dPCr</NU><DE>d<IT>t</IT></DE></FR> = <IT>F</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> + <FR><NU><IT>R</IT><SUB>PCr</SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (PCr<SUB>i</SUB> − PCr) (80)

<FR><NU>dCr</NU><DE>d<IT>t</IT></DE></FR> = <IT>F</IT><SUB>CK</SUB><IT>v</IT><SUB>CK</SUB> + <FR><NU><IT>R</IT><SUB>Cr</SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (Cr<SUB>i</SUB> − Cr) (81)

<FR><NU>dP<SUB>i</SUB></NU><DE>d<IT>t</IT></DE></FR> = <IT>F</IT><SUB>H</SUB><IT>H</IT><SUB>ATP</SUB>(<IT>t</IT>) + <FR><NU><IT>R</IT><SUB>P<SUB>i</SUB></SUB></NU><DE><IT>L</IT><SUB>m</SUB></DE></FR> (P<SUB>i,i</SUB> − P<SUB>i</SUB>) (82)
where Lm is the distance between the core of the myofibril and the mitochondrial membrane. In addition to Eqs. 77-82, the complete zero-dimensional model includes Eqs. 25-30 and 37-42.


    ACKNOWLEDGEMENTS

The authors thank Prof. Juri Engelbrecht (Tallinn, Estonia) for continuous support of this work, Drs. Bernard Korzeniewski (Krakow, Poland) and Maiys Aliev (Moscow, Russia) for very instructive discussions, and Dr. Bernard Korzeniewski for providing the computer program of his model. We have always greatly appreciated our most interesting and fruitful discussions with the late Prof. John R. Williamson (Johnson Foundation, Univ. of Pennsylvania, Philadelphia, PA), on the problems of modeling heart metabolism, for which he supplied the most reliable experimental data.


    FOOTNOTES

This work was supported by Estonian Science Foundation Grant 3204.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. §1734 solely to indicate this fact.

Address for reprint requests and other correspondence: V. Saks, Laboratory of Bioenergetics, Joseph Fourier University, BP 53X-38 041 Grenoble Cedex, France (E-mail: Valdur.Saks{at}ujf-grenoble.fr).

Received 21 July 1999; accepted in final form 10 November 1999.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
GLOSSARY
METHODS OF MODELING
RESULTS
DISCUSSION
REFERENCES
APPENDIX

1.   Aliev, MK, and Saks VA. Compartmentalized energy transfer in cardiomyocytes: use of mathematical modelling for analysis of in vivo regulation of respiration. Biophys J 73: 428-445, 1997[Abstract].

2.   Balaban, RS, Kantor HL, Katz LA, and Briggs RW. Relation between work and phosphate metabolite in the in vivo paced mammalian heart. Science 232: 1121-1123, 1986[ISI][Medline].

3.   Bereiter-Hahn, J, and Voth M. Dynamics of mitochondria in living cells: shape changes, dislocations, fusion and fission of mitochondria. Microsc Res Tech 27: 198-219, 1994[ISI][Medline].

4.   Brandes, R, and Bers D. Intracellular Ca increases the mitochondrial NADH concentration during elevated work in intact cardiac muscle. Circ Res 80: 82-87, 1997[Abstract/Free Full Text].

5.   Braunwald, E, Ross J, and Sonnenblick EH. Mechanisms of Contraction of the Normal and Failing Heart. Boston, MA: Little Brown, 1968.

6.   Brown, PN, Byrne GD, and Hindmarsh AC. VODE: a variable coefficient ODE Solver. SIAM J Sci Stat Comput 10: 1038-1051, 1989[ISI].

7.   Bruaset, AM, and Langtagen HP. A comprehensive set of tools for solving partial differential equations: Diffpack. In: Numerical Methods and Software Tools in Industrial Mathematics, edited by Daehlen M, and Tveito A.. Boston, MA: Birkhauser, 1997, p. 61-90.

8.   Cannell, MB, Cheng H, and Lederer WJ. The control of calcium release in heart muscle. Science 268: 1045-1049, 1995[ISI][Medline].

9.   Daut, J. The living cell as energy-transducing machine. A minimal model of myocardial energy metabolism. Biochim Biophys Acta 895: 41-62, 1987[ISI][Medline].

10.   Duszynski, J, Bogucka K, and Wojtczak L. Homeostasis of the protonmotive force in phosphorylating mitochondria. Biochim Biophys Acta 767: 540-547, 1984[ISI][Medline].

11.   Dzeja, P, Zeleznikar RJ, and Goldberg ND. Adenylate kinase: kinetic behavior in intact cells indicates it is integral to multiple cellular processes. Mol Cell Biochem 184: 169-182, 1998[ISI][Medline].

12.   Elzinga, G, and van der Laarse WJ. MVO2 max of the heart cannot be determined from uncoupled myocytes. Basic Res Cardiol 85: 315-317, 1990[ISI][Medline].

13.   Fossel, E, and Hoefeler H. A synthetic functional metabolic compartment. The role of propinquity in a linked pair of immobilised enzymes. Eur J Biochem 170: 165-171, 1987[Abstract].

14.   Hafner, RP, Brown GC, and Brand MD. Analysis of the control of respiration rate, phosphorylation rate, proton leak rate and protonmotive force in isolated mitochondria using the "top-down" approach of metabolic control theory. Eur J Biochem 188: 313-319, 1990[Abstract].

15.   Hamman, BL, Bittl JA, Jacobus WE, Allen PD, Spencer RS, Tian R, and Ingwall JS. Inhibition of creatine kinase reaction decreases the contractile reserve of isolated rat hearts. Am J Physiol Heart Circ Physiol 269: H1030-H1036, 1995[Abstract/Free Full Text].

16.   Hansford, RG, and Zorov D. Role of mitochondrial calcium transport in the control of substrate oxidation. Mol Cell Biochem 184: 359-369, 1998[ISI][Medline].

17.   Harrisson, GI, Wijhe MH, Groot de B, Dijk FJ, and van Beek JHGM CK inhibition accelerates transcytosolic energy signaling during rapid workload steps in isolated rabbit hearts. Am J Physiol Heart Circ Physiol 276: H134-H140, 1999[Abstract/Free Full Text].

18.   Hochachka, PW, and McClelland GB. Cellular metabolic homeostasis during large-scale change in ATP turnover rates in muscles. J Exp Biol 200: 381-386, 1997[Abstract/Free Full Text].

19.   Jacobus, WE, and Saks VA. Creatine kinase of heart mitochondria: changes in its kinetic properties induced by coupling to oxidative phosphorylation. Arch Biochem Biophys 219: 167-178, 1982[ISI][Medline].

20.   Kammermeier, H. Why do cells need phosphocreatine and a phosphocreatine shuttle? J Mol Cell Cardiol 19: 115-118, 1987[ISI][Medline].

21.   Korzeniewski, B. Regulation of ATP supply during muscle contraction: theoretical studies. Biochem J 330: 1189-1195, 1998[ISI][Medline].

22.   Korzeniewski, B, and Froncisz W. An extended model of oxidative phosphorylation. Biochim Biophys Acta 1060: 210-223, 1991[ISI][Medline].

23.   Meyer, RA, Sweeney HL, and Kushmerick MJ. A simple analysis of the "phosphocreatine shuttle." Am J Physiol Cell Physiol 246: C365-C377, 1984[Abstract].

24.   Mitchell, P. Compartmentation and communication in living systems. Ligand conduction: a general catalytic principle in chemical, osmotic and chemiosmotic reaction systems. Eur J Biochem 95: 1-20, 1979[Abstract].

25.   Mootha, VK, Arai AE, and Balaban RS. Maximum oxidative phosphorylation capacity of the mammalian heart. Am J Physiol Heart Circ Physiol 272: H769-H775, 1997[Abstract/Free Full Text].

26.   Nicolay, K, van der Toorn A, and Dijkhuizen RM. In vivo diffusion spectroscopy. An overview. NMR Biomed 8: 365-374, 1995[ISI][Medline].

27.   O'Rourke, B, Ramza BM, and Marban E. Oscillations of membrane current and excitability driven by metabolic oscillations in heart cells. Science 265: 962-966, 1994[ISI][Medline].

28.   Ovadi, J. Cell Architecture and Metabolic Channelling. New York: Springer-Verlag, 1995, p. 1-250.

29.   Rizzuto, R, Pinton P, Garrington W, Fay FS, Fogarty KE, Lifshitz LM, Tuft RA, and Pozzan T. Close contacts with the endoplasmic reticulum as determinants of mitochondrial Ca responses. Science 280: 1763-1766, 1998[Abstract/Free Full Text].

30.   Romashko, DN, Marban E, and O'Rourke B. Subcellular metabolic transients and mitochondrial redox waves in heart cells. Proc Natl Acad Sci USA 95: 1618-1623, 1998[Abstract/Free Full Text].

31.   Saks, VA, Chernousova GB, Gukovsky DE, Smirnov VN, and Chazov EI. Studies of energy transport in heart cells. Mitochondrial isoenzyme of creatine kinase: kinetic properties and regulatory action of Mg2+ ions. Eur J Biochem 57: 273-290, 1975[Abstract].

32.   Saks, VA, and Ventura-Clapier R (Editors). Cellular Bioenergetics: Role of Coupled Creatine Kinases. Boston, MA: Kluwer Academic, 1994, p. 1-346.

33.   Saupe, KW, Spindler M, Tian R, and Ingwall J. Impaired cardiac energetics in mice lacking muscle-specific isoenzymes of creatine kinase. Circ Res 82: 898-907, 1998[Abstract/Free Full Text].

34.   Veech, RL, Lawson JWR, Cornell NW, and Krebs HA. Cytosolic phosphorylation potential. J Biol Chem 254: 6538-6547, 1979[Abstract].

35.   Wan, B, Dounen C, Duszynsky J, Salama G, Vary TC, and Lanoue KF. Effect of cardiac work on electrical potential gradient across mitochondrial membrane in perfused hearts. Am J Physiol Heart Circ Physiol 265: H453-H460, 1993[Abstract/Free Full Text].

36.   Watchko, JF, Daood MJ, Sieck GC, LaBella JJ, Ameredes BT, Koretsky AP, and Wieringa B. Combined myofibrillar and mitochondrial creatine kinase deficiency impairs mouse diaphragm isotonic function. J Appl Physiol 82: 1416-1423, 1997[Abstract/Free Full Text].

37.   Wegmann, G, Zanolla E, Eppenberger HM, and Wallimann T. In situ compartmentation of creatine kinase in intact sarcomeric muscle: the acto-myosin overlap zone as a molecular sieve. J Muscle Res Cell Motil 13: 420-435, 1992[ISI][Medline].

38.   Williamson, JR, Ford G, Illingworth J, and Safer B. Coordination of citric acid activity with electron transport flux. Circ Res. 38, Suppl I: I39-I51, 1976[Medline].

39.   Zweier, JL, Jacobus WE, Korecky B, and Brandejs-Barry Y. Bioenergetic consequences of cardiac phosphocreatine depletion induced by creatine analog feeding. J Biol Chem 266: 20296-20304, 1991[Abstract/Free Full Text].


Am J Physiol Cell Physiol 278(4):C747-C764
0363-6143/00 $5.00 Copyright © 2000 the American Physiological Society