## Abstract

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 P_{i}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

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 (V˙o
_{2}) 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 ofV˙o
_{2} 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

### 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, P
_{i}) in myofibril and cytoplasm - Met
_{i} - Concentration of metabolite in mitochondrial IM space
- ATP
_{g}, ADP_{g} - ATP and ADP concentrations in microcompartment
- UQ
- Oxidized form of coenzyme Q
- c
^{3+} - Oxidized form of cytochrome
*c* - NAD
^{+} - NAD
^{+}in mitochondrial matrix - H
_{x} - H
^{+}in mitochondrial matrix - ATP
_{x} - ATP in mitochondrial matrix
- P
_{i,x} - P
_{i}in mitochondrial matrix

### Diffusion

- ∇
^{2} - Laplace operator
*D*_{Met}- Diffusion coefficient of metabolite
*G*_{CK}- Spatial distribution function of myofibrillar CK
*G*_{AK}- Spatial distribution function of myofibrillar AK
*G*_{H}- Spatial distribution function of myofibrillar ATPase

### Myofibrillar ATPase

*H*_{ATP}- Rate of ATP hydrolysis in myofibril
*H*_{ATP(max)}- Maximal
*H*_{ATP}

### CK in myofibril, myoplasm, and IM space

*v*_{CK}- CK reaction rate
*v*_{MiCK}*v*_{CK}in IM space*v*_{MiCK,G}*v*_{CK}, ATP and ADP from microcompartment*v*_{MiCK,I}*v*_{CK}, ATP and ADP from IM space*V*_{1},*V*_{−}_{1}- Maximal CK reaction rates in forward and reverse directions
*K*_{ia},*K*_{b},*K*_{ic},*K*_{d},*K*_{id},*K*_{ib},*K*_{Ib}- 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
*K*_{iaG},*K*_{icG}- MiCK dissociation constants for ATP and ADP in microcompartment

### AK in myofibrils and IM space

*v*_{AK}- AK reaction rate
*v*_{AKmit}*v*_{AK}in IM space*k*_{fa},*k*_{ba}- AK forward and backward reaction rate constants

### Mg^{2+}buffering

- mMet
- Mg
^{2+}-bound forms of metabolite (ATP or ADP) - fMet
- Mg
^{2+}-free forms of metabolite (ATP or ADP) - Mg
_{ext} - Mg
^{2+}in cytoplasm and myofibril - Mg
_{x} - Mg
^{2+}in mitochondrial matrix *K*_{DText},*K*_{DDext},*k*_{DTx},*k*_{DDx}- Mg
^{2+}dissociation constants

### H^{+} buffering

- H
_{ext}, H_{x} - H
^{+}concentration in cytoplasm and matrix - pH
_{ext}, pH_{x} - pH in cytoplasm and matrix
- Δ
*p*, ΔpH, Δψ - Protonmotive force, pH difference, and membrane potential across inner mitochondrial membrane
*r*_{buff},*r*_{buff,0}- H
^{+}buffering capacity coefficients

### Diffusion through mitochondrial outer membrane

**n⃗**- Unit vector, outward normal to the domain boundary
*L*_{i}- Width of the layer between inner and outer membranes
*L*_{m}- Width of the layer between core of myofibril and mitochondrial outer membrane
*R*_{Met}- Permeability of the mitochondiral outer membrane for the metabolite

### ATP diffusion between microcompartment and IM space

*v*_{exchATP},*v*_{exchADP}- Rate of ATP or ADP flux between microcompartment and IM space
*R*_{exchATP},*R*_{exchADP}- Rate constants

### ATP/ADP translocation by ANT

*v*_{ANT}- Net rate of ATP export by ANT from matrix
- Rate of ATP export by ANT from matrix to microcompartment
- Rate of ATP export by ANT from matrix to IM space
*V*_{ANT}- Maximal net rate of ATP export by ANT from matrix
*K*_{g},*K*_{i}- ANT reaction dissociation constants

### Fractional volumes

- FV
_{i} - Fractional volume of IM space with respect to cell volume
- FV
_{x} - Fractional volume of matrix with respect to cell volume

### Substrate dehydrogenation

*v*_{dh},*k*_{dh}- Net and maximum rates of substrate dehydrogenation
*K*_{mN}- Michaelis-Menten constant of substrate dehydrogenation for the NAD
^{+}/NADH ratio *p*_{D}- Relative sensitivity coefficient of substrate dehydrogenation to the NAD
^{+}/NADH ratio

### Phosphate carrier

*v*_{PI},*v*_{f},_{PI},*v*_{b},_{PI}- Net, forward (to matrix), and backward rates of phosphate carrier
- p
*K*_{a} - One-half dissociation constant for monovalent P
_{i}

### Proton leak

*k*_{L1},*k*_{L2}- Phenomenological constants for proton leak

### ATP synthesis

*v*_{sn},*k*_{sn}- Net and maximal rates of ATP synthase
- Δ
*G*_{0}, Δ*G*_{sn} - Gibbs free energy and thermodynamic span of ATP synthase
*n*_{a}- H
^{+}/ATP stoichiometry for ATP synthase

### Respiratory chain

*v*_{C1},*v*_{C3},- Rates of complexes I, III, and IV
*k*_{C1},*k*_{C3},*k*_{C4}- Maximum rates of complexes I, III, and IV
- a
^{2+}, c^{2+}, a^{3+}, c^{3+} - Concentrations of reduced and oxidized forms of cytochromes
*a*_{3}and*c* *k*_{MO}- Michealis-Menten constant of complex IV for oxygen
*E*_{n},*E*_{u},*E*_{c},*E*_{a}- Redox potentials of NAD, ubiquinone, and cytochromes
*c*and*a*_{3} *E*_{n,0},*E*_{u,0},*E*_{c,0},*E*_{a,0}- Standard redox potentials of NAD, ubiquinone, and cytochromes
*c*and*a*_{3} - A
_{x}, NAD_{tot}, O_{2}UQ_{tot},*c*_{tot},*a*_{tot} - Total concentrations of adenine nucleotides, NAD, oxygen in matrix, ubiquinone, and cytochromes
*c*and*a*_{3}. *r*_{buff},*r*_{buff,0}- Buffering capacity coefficients

### General constants

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

## METHODS OF MODELING

#### 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 P_{i}. 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 P_{i}. 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. P_{i} is transported by a phosphate carrier.

#### 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. 1
*B*. 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 V˙o
_{2}. 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 O_{2} ⋅ 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 O_{2}/ ⋅ 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 maximalV˙o
_{2}, 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 minimalV˙o
_{2} is 12–15 μmol O_{2} ⋅ 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 10^{5} 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 10^{5}-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
). 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 . All values of the experimentally determined parameters used in the model were taken from literature and are shown in Table 1.

#### 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

#### 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.

#### 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,V˙o
_{2} 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 V˙o
_{2} 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 maximalV˙o
_{2} recorded by Williamson et al. (38) is close to 160 μmol O_{2} ⋅ 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 whenV˙o
_{2} exceeds 80–100 μmol O_{2} ⋅ 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).

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, ofV˙o
_{2} (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 whenV˙o
_{2} 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 betweenV˙o
_{2} 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 P_{i} (Fig.9).

#### 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 correspondingV˙o
_{2} 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 maximalV˙o
_{2} 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 calculatedV˙o
_{2} 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 O_{2} ⋅ g dry wt^{−}
^{1} ⋅ min^{−}
^{1}by using the values of the diffusion coefficients characteristic of water. The dotted line in Fig. 10
*A* was obtained for the case of infinitely fast diffusion of ADP (zero-dimensional model, seemethods 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. 10
*A*). 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.10
*A*).

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. 10
*B*). 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 11
*A* 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.10
*A* correspond to maximal and minimal ADP concentrations in Fig. 11
*A*. Figure 11
*B* 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. 10
*B*, in this case, one observes very remarkable ADP concentration gradients within the cell.

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 P_{i} 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 P_{i} 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.

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 P_{i}; 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%.

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

The results of this study show that the linear dependence ofV˙o
_{2} 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 ofV˙o
_{2} 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 and11). 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 V˙o
_{2}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 ofV˙o
_{2} 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 P_{i}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 V˙o
_{2} 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 P_{i} 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.

## Acknowledgments

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.

## 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
Equation 1
Equation 2
Equation 3
Equation 4
Equation 5
Equation 6 where ATP, ADP, AMP, PCr, Cr, and P_{i} are the total concentrations of the metabolites in the myofibril and myoplasm,*D*
_{Met} is the diffusion coefficient of the metabolite (Met, i.e., ATP, ADP, AMP, PCr, Cr, P_{i}),*v*
_{CK} is the CK reaction rate,*v*
_{AK} is the AK reaction rate, and*H*
_{ATP} is the rate of ATP hydrolysis in the myofibril. The factors *G*
_{CK},*G*
_{AK}, and *G*
_{H} describe the spatial inhomogeneities of myofibrillar CK, AK, and ATPase, respectively. The symbol ∇^{2} denotes the Laplace operator.

The *v*
_{CK} is described by *Eq. 7
*
Equation 7 where
Equation 8and
Equation 9 where mATP and mADP are the concentrations of the Mg^{2+}-bound forms of ATP and ADP, respectively, and fATP and fADP are the Mg^{2+}-free forms. The total concentrations of adenylates and their Mg^{2+}-bound forms are related as follows
Equation 10
Equation 11
Equation 12
Equation 13 where Mg is the level of free Mg^{2+} and *K*
_{DT}and *K*
_{DD} are Mg^{2+} dissociation constants. The level of free Mg^{2+} and its dissociation constant for the matrix differ from their corresponding values for the rest of the cell (22). The *v*
_{AK} is described as follows
Equation 14 where*k*
_{fa} and *k*
_{ba} are the forward and backward reaction rate constants for the AK reaction.

We approximate the hydrolysis rate by a piecewise linear time-periodic function *H*
_{ATP}. The rate increases from 0 to*H*
_{ATP(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
Equation 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
Equation 16the distribution of AK by
Equation 17 and the distribution of ATPase (5) by
Equation 18where *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*G*
_{CK} and *G*
_{AK} over the myofibril and myoplasm (0 ≤ *x* ≤ 1, 0 ≤ *y *≤ 1.2) and*G*
_{H} 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
Equation 19
Equation 20
Equation 21
Equation 22
Equation 23
Equation 24 where Met_{i} is the concentration of the metabolite in the mitochondrial IM space and *R*
_{Met} 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 P_{i} transportation through the inner membrane according to the following equations
Equation 25
Equation 26
Equation 27
Equation 28
Equation 29
Equation 30 where*v*
_{MiCK} is MiCK reaction rate,*v*
_{ANT} is net rate of ATP export by ANT from the matrix, *v*
_{AKmit} is AK reaction rate (see *Eq.14
*), *v*
is rate of the phosphate carrier in the inner membrane,*L*
_{i} is width of the layer between the inner and outer membranes, and FV_{i} is fractional volume of the IM space with respect to the cell volume.

The MiCK reaction rate (*v*
_{MiCK}) is divided into*v*
_{MiCK,G}, which involves microcompartment ATP(ADP), and *v*
_{MiCK,I}, which involves ATP(ADP), from the IM space
Equation 31 The reaction rates are as follows
Equation 32
Equation 33where

Equation 34

Restricted diffusion between the microcompartment and the IM space is governed by Equation 35 for ATP and by Equation 36 for ADP.

For the oxidative phosphorylation, we use the model by Korzeniewski (21)
Equation 37
Equation 38
Equation 39
Equation 40
Equation 41
Equation 42 where UQ is the oxidized form of coenzyme Q, c^{3+} is the oxidized form of cytochrome *c,* NAD^{+} is NAD^{+} in the matrix, H_{x}, ATP_{x}, and P_{i,x} are H^{+}, ATP, and P_{i} in the matrix, and FV_{x} is fractional volume of the mitochondrial matrix with respect to the cell volume. The following reaction rates are used in the above equations (21)
Equation 43for substrate dehydrogenation (dh)
Equation 44 for complex I (C1)
Equation 45for complex III (C3)
Equation 46for complex IV (O_{2})
Equation 47 for ATP synthase (sn)
Equation 48 for phosphate carrier (PI), and
Equation 49for proton leak (leak). In *Eqs. 43-49
*, we used the total (tot) metabolite pools for substrate (NAD), coenzyme Q (UQ), cytochrome *c*(c), cytochrome *a*
_{3} (a), and adenine nucleotides (in matrix, A_{x}) as follows
Equation 50
Equation 51
Equation 52
Equation 53
Equation 54 The concentrations of Mg^{2+}-bound and Mg^{2+}-free forms of adenine nucleotides in the matrix are obtained using *Eqs.10-13
*. For computing pH in the matrix (pH_{x}), pH difference across the inner membrane (ΔpH), protonmotive force (Δ*p*), electrical potential (Δψ), and buffering rate (*r*
_{buff}), the following equations were used
Equation 55
Equation 56
Equation 57
Equation 58
Equation 59
Equation 60 where the coefficient *u* = ΔΨ/Δ*p* is kept constant. Because of the high permeability of the mitochondrial outer membrane for protons, pH in IM space (pH_{ext}) is considered to be equal to that in the cytoplasm and is considered to be constant. The coefficient Ζ is used to convert between molar concentrations and the concentrations used in the model (μM). The coeffiient *Z* is expressed as
Equation 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
Equation 62 where*v*
is the rate of ATP export by ANT from the matrix to the microcompartment and*v*
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
Equation 63
Equation 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
Equation 65
Equation 66 *Equations65
* and *
66
* can be solved for ATP_{g} and ADP_{g} by using, for example, the Newton-Raphson method.

The thermodynamic span of ATP synthase is calculated as follows
Equation 67 where*n*
_{a} and Δ*G*
_{0} 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
Equation 68
Equation 69
Equation 70 respectively. For computing the concentration of reduced cytochrome*a*
_{3} used in *Eq. 46
*, the redox potential of cytochrome *a*
_{3}
Equation 71 is substituted into the relation
Equation 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
Equation 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 *G*
_{CK},*G*
_{AK}, and *G*
_{H} leads us to
Equation 74
Equation 75 and
Equation 76

#### Zero-dimensional model.

The concentrations of the free metabolites in cardiac myocytes are described by the following equations
Equation 77
Equation 78
Equation 79
Equation 80
Equation 81
Equation 82 where*L*
_{m} 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
*.

## Footnotes

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).

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.

- Copyright © 2000 the American Physiological Society