Introduction

The energy demand imposed by physical exercise on the components of the oxygen transport and utilization system requires a close link between cellular and external respiration in order to maintain adenosine triphosphate (ATP) homeostasis. At the onset of exercise, an immediate increase in the rate of production of (ATP) in active skeletal muscle fibers is required to meet the increased rate of ATP utilization (metabolic demand). Mechanical power output of muscle – measured on a cycle ergometer – responds quickly (τPower < 2 s to a step change in work rate, indicating an even faster response of the ATP production rate by contracting muscles. In contrast, the corresponding dynamic responses of pulmonary oxygen uptake (VO2p) and muscle oxygen uptake (VO2m)29,38,47,50,55 are much slower \({(\tau_{{\rm VO}_{2}} \sim 30\,\hbox{s})}\). This slower response reflects the ATP production rate from oxidation of reducing equivalents in active muscle. While the steady states for ATP demand and ATP supply through oxidative metabolism match perfectly after ∼2–3 min at submaximal work rates, their transient responses differ by at least an order of magnitude.

Noninvasive measurements of VO2p dynamics can characterize differences in onset kinetics of healthy and chronically ill individuals.7,8,49 These measurements, however, not only reflect cellular metabolism, but also the effect of components of the delivery system that transports oxygen from the external environment to the mitochondria (e.g., pulmonary ventilation, diffusion from alveoli to blood, transport in blood – dissolved and bound to hemoglobin – to tissues, diffusion from blood to cells). Because measurements of VO2p dynamics do not accurately indicate the cellular O2 availability and utilization rate in contracting muscle, they are insufficient for deducing mechanisms of metabolic control at the cellular level.

To investigate factors controlling the rate of pulmonary oxygen uptake at exercise onset, oxygen uptake has been measured directly across the femoral bed during cycling in humans29 and across isolated whole skeletal muscle preparations in dogs.24–27,29 These whole-muscle studies have provided valuable insights about the effect of altering convective oxygen delivery and peripheral oxygen diffusion on the dynamics of muscle oxygen uptake during moderate and heavy intensity exercise. Nevertheless, they also have some inherent limitations. When a whole muscle is stimulated in lieu of voluntary muscle contraction, the fiber-type recruitment pattern is different. Furthermore, controlling bulk oxygen delivery does not necessarily control the matching of perfusion to metabolism. Measurements of the oxygen content in the arterial and venous blood do not reflect the actual amount of oxygen available at the mitochondrial level. Local oxygen availability (i.e., cellular PO2) in intact muscle tissue is difficult to assess during exercise in these whole-tissue preparations.

Cellular oxygen availability can be assessed in muscle of contrasting fiber types with phosphorescence quenching9 or in intact working muscle with either near-infrared spectroscopy28 or 1H magnetic resonance spectroscopy.36,60 However, these methods lack specificity when measuring the dynamics of cellular oxygen consumption in a large muscle group during exercise. Nevertheless, intracellular oxygen consumption can be estimated using 31P magnetic resonance spectroscopy by using PCr dynamics as a proxy (within 10%) to muscle oxygen consumption dynamics.3,53,54,65 A more suitable approach to investigate mechanisms controlling respiration is to conduct simultaneous measurements of putative intramuscular control mediators (PCr, Cr, Pi, ATP, ADP) dynamics and VO2p dynamics while performing exercise involving a large muscle group (e.g., quadriceps). This approach provides a large change in VO2p at exercise onset and reliable VO2p dynamic responses which can be compared to those of PCr.53,65

Even if measurements at each scale (i.e., organelle, cellular, tissue/organ, and whole organism level) were experimentally feasible with the appropriate time resolution, these must be integrated into a coordinated system for data interpretation and for elucidation of mechanisms of metabolic control and regulation. For quantitative understanding, integration of transport and metabolic processes at the cellular, tissue/organ, and whole-body levels requires a formal theoretical framework. As a complement to experimental studies of oxygen uptake/consumption at various exercise intensities and biological levels, we used a mathematical model to integrate data and to simulate responses at the cellular and tissue level. The models integrated in this study have been validated previously.42,62

The “systems biology” approach employed here in linking external to internal respiration uses both experimental data from non-invasive methods and computational models to understand complex phenomena and biological systems underlying physiological function and ATP homeostasis. For initial model development to analyze exercise responses, we used data at the whole body and tissue level. In addition, we measured pulmonary oxygen uptake (indirect calorimetry) and muscle oxygenation (near-infrared spectroscopy) dynamics during exercise in healthy volunteers.42 Based on these data, we developed a computational model of oxygen transport and metabolism in skeletal muscle that can simulate and predict muscle oxygen consumption dynamics during exercise.42 In this study, we integrated our model with a model of energy balance in myocytes41,62 to provide a more mechanistic expression for oxidative phosphorylation,35 which is linked to the ATPase and creatine kinase (CK) reactions.62 With this model, muscle oxygen consumption dynamics at the cellular level can be estimated during exercise. Finally, model simulations are used to provide insight with respect to concentration and flux rate dynamics of key metabolites (PCr, ATP, ADP, O2) participating in oxidative phosphorylation and ATP homeostasis.

Methods

Model Development

Oxygen utilization rate in skeletal muscle (UO2m) is linked to oxygen uptake rate in lungs (VO2p) by transport processes within the tissues and via the cardiovascular system (Fig. 1). To analyze the dynamic relationship between UO2m and VO2p, a mathematical model is needed. In an earlier approach to developing this relationship, simulations were based on a model42 that approximated oxygen utilization in working skeletal muscle as a mono-exponential function of work rate.10 This model was used to analyze dynamic responses to exercise at several levels of intensity.42 With this simple model, however, distinct metabolic processes involved during exercise stimulus cannot be evaluated. Therefore, a more mechanistic model is needed that describes oxidative phosphorylation by incorporating ATPase and CK reactions. Such a metabolic model has been developed previously.62 Here, this model of oxidative phosphorylation in myocytes is combined with a model of oxygen transport and utilization42 for analysis of responses to exercise.

Figure 1
figure 1

Oxygen utilization and transport between lungs and skeletal muscle and cellular metabolism during exercise

Full size image

Mass Transport Balances

Oxygen concentration dynamics in skeletal muscle are represented by compartmental mass balances. These balance equations involve total oxygen (C Tcap , C Ttis ) to free oxygen (C Fcap , C Ftis ) in the capillaries and tissues:

$$ V_{\rm cap} \frac{dC_{\rm cap}^{\rm T}}{dt}=Q_{\rm m} (t)\left( C_{\rm art}^{\rm T} -C_{\rm cap}^{\rm T} \right)-\hbox{PS}(t)\left( C_{\rm cap}^{\rm F} -C_{\rm tis}^{\rm F} \right) $$
(1)
$$ V_{\rm tis} \frac{dC_{\rm tis}^{\rm T}}{dt}=\hbox{PS}(t)\left(C_{\rm cap}^{\rm F} -C_{\rm tis}^{\rm F}\right)-\hbox{UO}_{\rm 2m} (t) $$
(2)

Based on relationships between total, bound, and free oxygen concentrations, these equations have an equivalent representation with derivatives of free oxygen (Appendix I). The muscle blood flow, Q m(t), the capillary-tissue transport coefficient represented as the permeability – surface area product, PS(t), and the oxygen utilization, UO2m(t), depend on exercise intensity. The rate of muscle oxygen uptake is defined as \({\hbox{VO}_{\rm 2m} (t)=Q_{\rm m} (t)\left( C_{\rm art}^{\rm T} -C_{\rm cap}^{\rm T} \right)}\) (All variables and symbols are defined in the “List of Symbols”).

Metabolic Reaction Balances

The metabolic reaction processes of oxidative phosphorylation during an exercise stimulus have been reported in Fig. 1 and can be expressed as the concentration dynamics of ATP and PCr:

$$ \frac{dC_{\rm ATP}}{dt}=-\phi_{\rm ATPase} +\beta \phi_{\rm OxPhos} -\phi_{\rm CK}^{\rm r}+\phi_{\rm CK}^{\rm f} $$
(3)
$$ \frac{dC_{\rm PCr}}{dt}=\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f} $$
(4)

where β is the P:O2 ratio in vivo and the reaction fluxes φ j are functions of the ATP and PCr concentrations, which must satisfy the mass balances for conservation of adenosine and creatine:

$$ C_{\rm A,tot} =C_{\rm ADP} +C_{\rm ATP} $$
(5)
$$ C_{\rm Cr,tot} =C_{\rm Cr} +C_{\rm PCr} $$
(6)

The metabolic flux of oxygen is nonlinear related to the ADP and oxygen concentrations:

$$ \phi_{\rm OxPhos} =V_{\rm max}^{\rm j} \left( \frac{C_{\rm ADP}}{K_{\rm ADP} +C_{\rm ADP}} \right)\left( \frac{C_{\rm tis}^{\rm F}}{K_{\rm m} +C_{\rm tis}^{\rm F}} \right) $$
(7)

where the rate coefficient, V jmax , depends on exercise intensity. The metabolic flux for the ATPase reaction is proportional to the ATP concentration:

$$ \phi_{\rm ATPase} =k_{\rm ATPase}^{\rm j} C_{\rm ATP} $$
(8)

where the reaction rate coefficient, k jATPase , depends on exercise intensity. The forward and reverse reaction fluxes of CK are nonlinearly related to the coupled concentrations of Cr, PCr, ADP, and ATP:62

$$ \phi_{\rm CK}^{\rm r} =\frac{V_{\rm CK}^{\rm r} ({\frac{{C_{\rm Cr} C_{\rm ATP}}}{{K_{\rm iq}K_{\rm p}}}})} {1+\frac{C_{\rm ADP}}{K_{\rm ia}}+\frac{C_{\rm ATP}}{K_{\rm iq}} +\frac{C_{\rm PCr}}{K_{\rm ib}}+\frac{C_{\rm ADP} C_{\rm PCr}}{K_{\rm b} K_{\rm ia}} +\frac{C_{\rm Cr} C_{\rm ATP}}{K_{\rm iq} K_{\rm p}}} $$
(9)
$$ \phi_{\rm CK}^{\rm f} =\frac{V_{\rm CK}^{\rm f} \frac{C_{\rm ADP} C_{\rm PCr}}{K_{\rm b} K_{\rm ia}}} {1+\frac{C_{\rm ADP}}{K_{\rm ia}}+\frac{C_{\rm ATP}}{K_{\rm iq}} +\frac{C_{\rm PCr}}{K_{\rm ib}}+\frac{C_{\rm ADP} C_{\rm PCr}} {K_{\rm b} K_{\rm ia} }+\frac{C_{\rm Cr} C_{\rm ATP}}{K_{\rm iq} K_{\rm p}}} $$
(10)

The contributions of glycogenolysis and glycolysis to ATP synthesis are not included in the mathematical model Eq. (3) at this stage of model development. Thus, results obtained from simulations with this incomplete model need to be interpreted accordingly.

The reaction processes, Eqs. (3) and (4), are related to the transport processes, Eqs. (1) and (2), through the rates of oxygen metabolism, φOxPhos, and oxygen utilization:

$$ \phi_{\rm OxPhos} V_{\rm tis}=\hbox{UO}_{\rm 2m}^{\rm j} $$
(11)

which involves a tissue volume, V tis, and depends on the exercise intensity j. The total muscle volume is defined as \({V_{\rm m} =V_{\rm cap} +V_{\rm tis}}\) where V m is computed as a fraction of the whole body mass (49%) and V cap, V tis are each computed as a fraction of the V m.42

Exercise-dependent Functions

In response to a step increase in work rate from a steady-state warm up condition, the dynamic response of blood flow Q jm at exercise intensity j is assumed to be exponential:42

$$ Q_{\rm m}^{\rm j} (t)=Q_{\rm m}^{\rm W} +\Delta Q_{\rm m}^{\rm j} \left[ 1-\hbox{exp} \left( t^{\rm W}-t \right)/\tau _{Q_{\rm m}}^{\rm j} \right] $$
(12)

where Q Wm is the steady-state value during warm-up, ΔQ jm is the increase in blood flow, and \({\tau_{Q_{\rm m}}^{\rm j}}\) is the time constants of muscle blood flow, and t W is the initial time. In response to exercise, blood flow increases, which increases the rate of capillary-tissue transport:11

$$ \hbox{PS}^{\rm j}(t)=\hbox{PS}^{\rm R}+\Delta \hbox{PS} \left[ 1-\hbox{exp} \left( Q_{\rm m}^{\rm R} -Q_{\rm m}^{\rm j} \right)/Q_{\rm C} \right] $$
(13)

where \({\Delta \hbox{PS}=\hbox{PS}^{\rm E}-\hbox{PS}^{\rm R}}\), PSR and Q C are constants. For a description of their meaning refer to the List of Symbol.

Model Simulation

For comparison of simulated and experimental responses of oxygen obtained from human exercise experiments at several intensities, muscle oxygen saturation is computed:42

$$ \hbox{StO}_{\rm 2m} =\left[ C_{\rm cap}^{\rm B} V_{\rm cap} +C_{\rm tis}^{\rm B} V_{\rm tis} \right]/T_{\rm m} $$
(14)

The parameters T m (total amount of Hb and Mb), V cap and V tis are constants independent of exercise intensity. The initial conditions for the simulations of the dynamic responses are specified at a warm-up steady state (Appendix II):

$$ t=t^{\rm W},\quad C_{\rm cap}^{\rm F} =C_{\rm cap}^{\rm F,W},\quad C_{\rm tis}^{\rm F} =C_{\rm tis}^{\rm F,W},\quad C_{\rm ATP} =C_{\rm ATP}^{\rm W},\quad C_{\rm PCr} =C_{\rm PCr}^{\rm W} $$
(15)

To simulate responses to exercise, the model equations were solved numerically using a robust algorithm for stiff ordinary differential equations (DLSODE, http://www.netlib.org/odepack/).31

Parameter Estimation

Values of most parameters in the model equations above have been determined previously and reported in Tables 1A, B and 2.42,62 Others are evaluated directly from steady-state relationships (k jATPase , PS R) (Table 3) or by optimal estimation (ΔQ jm , V jmax ) (Table 4).

Table 1 (A) Values of muscle oxygen utilization, UO j2m , and pulmonary oxygen uptake, VO j2p , (mmol min−1) at various exercise intensities (j = R,W,M,H,V)42 and (B) Values of blood flow Q jm (L min−1) at various exercise intensities (j = R,W), mean response times τ jQ_m (s) (j = M,H,V) and compartment volumes V(L) for model simulations.42
Full size table
Table 2 Model parameters values for simulation of all experiments.42,62
Full size table
Table 3 Means (n = 7) of reaction rate coefficients for ATPase k jATPase (min−1) at different exercise intensities (j = R,W,M,H,V) and of permeability-surface area at rest (L min−1)
Full size table
Table 4 Mean (n = 7) parameter values for exercise intensity j = M,H,V: change in muscle blood flow ΔQ jm (L min−1) from warm-up; maximal rate of oxidative phosphorylation V jmax (mM min−1)
Full size table

From experiments, the values of V tis and UO j2m are determined so that at each exercise intensity we can evaluate φOxPhos from Eq. (11). Under steady-state conditions and assuming β = 6 (32), Eqs. (3) and (4) simplify to

$$0=-\phi_{\rm ATPase}+6\phi_{\rm OxPhos} $$
(16)

Combining this with Eq. (8) and given a constant value of C ATP, due to the ATP homeostasis, we can evaluate k jATPase (see Table 3). To calculate PSR, we consider Eqs. (1) and (2) at steady state under resting conditions:

$$ 0=Q_{\rm m}^{\rm R} \left( C_{\rm art}^{\rm T} -C_{\rm cap}^{\rm T} \right)-\hbox{PS}^{\rm R}\left( C_{\rm cap}^{\rm F}-C_{\rm tis}^{\rm F} \right) $$
(17)
$$ 0=\hbox{PS}^{\rm R}\left( C_{\rm cap}^{\rm F} -C_{\rm tis}^{\rm F} \right)-\hbox{UO}_{\rm 2m}^{\rm R} $$
(18)

From experiments, C Tart , C F,Rtis , Q jm and UO R2m are determined. By simultaneous solution of Eqs. (17) and (18), incorporating the relationship \({C_{\rm cap}^{\rm T,R}(C_{\rm cap}^{\rm F,R})}\), we estimate PSR and C F,Rcap .

The parameters ΔQ jm and V jmax must be estimated for each subject at each exercise intensity. For this purpose, we find the parameter values that yield the best fit of the model output dynamics (from numerical solution of the model equations) to the experimental data. Specifically, for each subject at each exercise intensity, we minimize a least-squares objective function:

$$ \Gamma (\Delta Q_{\rm m}^{\rm j} ,V_{\rm max}^{\rm j})=\frac{1}{2}\sum\limits_{i=1}^N \left[ \left( \frac{\hbox{StO}_{\rm 2m} (t_i )-\hbox{StO}_{\rm 2m}^{\rm W}}{\hbox{StO}_{\rm 2m}^{\rm W}} \right)_{\rm exp} -\left( \frac{\hbox{StO}_{\rm 2m} (t_i )-\hbox{StO}_{\rm 2m}^{\rm W}}{\hbox{StO}_{\rm 2m}^{\rm W}} \right)_{\rm mod} \right]^2 $$
(19)

where N is the number of data points. The objective function is minimized by numerical optimization using adaptive, non-linear algorithm (DN2FB, http://www.netlib.org).19

Results

In response to different exercise intensities, simulations represent dynamic responses of muscle ATP and PCr concentrations, oxygen saturation, StO2m, ATPase (φATPase), oxidative phosphorylation (φOxPhos), and net CK \({(\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f})}\) flux rates. By fitting of simulated StO2m to experimental data, optimal estimates were obtained for the maximal flux rate of oxidative phosphorylation V jmax and change in muscle blood flow ΔQ jm of seven subjects at exercise intensities j = M,H,V. From model simulations, mean response times \({(\tau_{{\rm VO}_{\rm 2m}}, \tau_{{\rm UO}_{\rm 2m}}, \tau_{\rm PCr})}\) were computed for muscle O2 uptake, muscle O2 utilization, and PCr, respectively.

Typical StO2m step responses to a change in exercise intensity relative to a warm-up steady state, StO W2m , for one subject are presented in Fig. 2. The decrease in StO2m is greater for a higher exercise intensity. Referring to the same subject, model simulations of UO j2m responses to increased exercise intensities (j = M,H,V) are compared to measured \({\hbox{VO}_{\rm 2p}^{\rm j}}\) (Fig. 3). Dynamic responses of \({\hbox{UO}_{\rm 2m}^{\rm j}}\) are faster than the dynamic responses of \({\hbox{VO}_{\rm 2m}^{\rm j}}\) for all subjects. The corresponding mean response times are compared in Table 5.

Figure 2
figure 2

Relative oxygen saturation in muscle, StO j2m : for representative subject responses to step changes from a warm-up steady-state condition (W) to a steady state during moderate, heavy, and very heavy exercise (j = M,H,V). Model output compared with experimental data.42

Full size image
Figure 3
figure 3

Comparison between dynamic responses of experimental pulmonary oxygen uptake and simulated muscle oxygen consumption for a representative subject.42 Forcing functions are step changes from a warm-up steady-state condition (W) to a work rate of moderate, heavy, and very heavy intensity exercise (j = M,H,V). Muscle oxygen consumption and pulmonary oxygen uptake are represented (a) in absolute terms as VO j2p and UO j2m , respectively and (b) as absolute changes from warm-up steady-state values, i.e., \({\hbox{VO}_{\rm 2p}^{\rm j} -\hbox{VO}_{\rm 2p}^{\rm W}}\) and \({\hbox{UO}_{\rm 2m}^{\rm j} -\hbox{UO}_{\rm 2m}^{\rm W}}\), respectively

Full size image
Table 5 Mean (N = 7) response times for exercise intensity j = M,H,V of muscle oxygen uptake \({{\varvec \tau}_{{\rm VO}_{\rm 2m}}^{\rm j}}\) (s), muscle oxygen utilization \({{\varvec \tau}_{{\rm UO}_{\rm 2m}}^{\rm j}}\) (s) and PCr dynamics τ jPCr (s)
Full size table

From optimal estimation of V jmax from seven subjects, their mean values are independent of exercise intensity (Table 4). In contrast, the mean change in muscle blood flow ΔQ jm increased by about 80% with higher exercise intensity.

Model simulated metabolic flux rates φOxPhos and \({\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f} }\) in response to a step change in φATPase (Table 3), corresponding to a given exercise intensity, are presented in Fig. 4. As expected (Eq. 3), at all exercise intensities, 6φOxPhos reaches the same steady as φ ATPase. Accordingly to the direct proportionality of φOxPhos and \({\hbox{UO}_{\rm 2m}^{\rm j}}\) (Eq. 11), the mean response time of φOxPhos given by \({\tau_{{\rm UO}_{\rm 2m}}}\) (Table 5) increases by about 35% from moderate to very high intensity exercise. While φOxPhos increases monotonically to the steady state, \({\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f}}\) shows an initial step decrease and then increases in exponential manner towards zero. The dynamic responses of the \({\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f}}\) and φOxPhos fluxes are closely balanced to guarantee ATP homeostasis.

Figure 4
figure 4

Simulation results for representative subject for the fluxes of ATPase, φ jATPase , oxygen phosphorylation, φ jOxPhos and creatine kinase,\({{\varvec{\phi}}_{\rm CK}^{\rm r,j}-{\varvec{\phi}}_{\rm CK}^{\rm f,j}}\). The simulations represent the dynamic between the warm-up steady state and the three levels of exercise intensity: moderate, heavy, and very heavy (j = M,H,V)

Full size image

The C ADP at steady-state, during exercise at all intensities, reached in our simulations (i.e., Figs. 2–5), consistently increases from moderate to very heavy intensity of exercise (M, 2.32 · 10−3 mM; H, 5.26 · 10−3 mM; V, 1.45 · 10−2 mM).

Figure 5
figure 5

Simulation results for representative subject. Comparison of absolute variation between the PCr break-down and the ATP consumption. The simulations are carried out starting from the warm-up steady state to the three levels of exercise intensity: moderate (continuous), heavy (dashed) and very heavy (dotted)

Full size image

The ATP concentration change is small compare to the PCr concentration change which decreases in response to a change in exercise intensity (Fig. 5). The mean response time of PCr, τPCr, shows a change with exercise that is similar to the change of \({\tau_{{\rm UO}_{\rm 2m}}}\). Model simulations were also performed assuming that the Michaelis-Menten parameter K m is sufficiently small in Eq. (7) such that φOxPhos is independent of oxygen concentration. In response to exercise, the mean response times for \({\hbox{UO}_{\rm 2m}^{\rm j}}\) and \({C_{\rm PCr} (\tau_{{\rm UO}_{\rm 2m}}, \tau_{\rm PCr})}\) are nearly the same for moderate, heavy and very heavy exercise (13, 13, 15 s). These response times are about 2–3 s less than those for which φOxPhos depends on oxygen concentration (Table 5). Also, for φOxPhos with no oxygen dependence, the mean values of V jmax (44 (M), 43 (H), 37 (V) mM min−1) compared to those given in Table 4 differ only for the very heavy intensity (j = V). Although V Vmax decreased from 44 to 37 mM min−1 without oxygen dependence, the mean response time of muscle oxygen uptake \({\tau_{{\rm VO}_{\rm 2m}}}\) did not change.

A reduction of the muscle volume, V m, from 40 to 9 l leads to greater amplitude changes of φATPase, net \({\phi_{\rm CK}^{\rm r}-\phi_{\rm CK}^{\rm f}}\) fluxes, C PCr and about a 50% increase in τPCr and \({\tau_{{\rm UO}_{\rm 2m}}}\). Also, the estimated V jmax in this range lies between 45 and 25 mM min−1 and was independent of exercise intensity.

Discussion

A mathematical model based on mass balances and metabolic and transport fluxes was successfully developed and applied to study the interacting regulation of cellular respiration and bioenergetics required to maintain ATP homeostasis during step transitions from rest to work rate intensities representing moderate, heavy, and very heavy exercise on a cycle ergometer in humans. This model couples oxygen transport and utilization in skeletal muscle42 to a model of cellular metabolism and energetics.62 The latter model, which describes oxidative phosphorylation, ATPase and CK reactions in finger and wrist flexors during an imposed stimulus, was adapted to investigate cellular energetics in large muscle groups.

In this study, transport and metabolic processes at the onset of exercise were quantified in the steady state and during the exercise transition by comparing model simulations with measurements of pulmonary oxygen uptake and muscle oxygen saturation at different exercise intensities on a cycle ergometer. Then, the dynamics of the metabolic flux rates of φOxPhos, φATPase, φ fCK and φ rCK were simulated for each of the seven subjects at three levels of exercise intensity.

Internal and External Respiration

The regulation of cellular oxygen consumption in skeletal muscle during exercise in vivo depends on electron flow, proton pumping, metabolic fluxes of NADH, ADP-dependent feedback control and oxygen delivery.16 At the cellular level, the dynamics of oxidative phosphorylation to a step increase in energy demand is expected to be fast (<450 ms) based on studies with isolated mitochondria in vitro.12,58 This response time is at least two orders of magnitude faster than in vivo measurements64,65 of oxygen uptake dynamics at the mouth (∼45 s). Since direct in vivo measurements of muscle oxygen consumption dynamics (UO2m) are difficult to obtain in humans during exercise, an indirect estimate is obtained from pulmonary oxygen uptake (VO2p) based on measurable variables at the airway opening.66 Under non-steady state conditions, however, VO2p cannot provide a reliable estimate of muscle oxygen consumption because of differences in dynamics at exercise onset.37 Experimentally, Grassi et al.29 found no significant differences in response dynamics between muscle and pulmonary uptake, VO2m and VO2p, during the transition from light to moderate intensity exercise. The methodology used for measuring VO2m dynamics across the femoral bed during exercise may not be sufficiently accurate because the VO2m response depends on lumped measurements of oxygen content in the femoral vein and bulk measurements of femoral blood flow,61 which may result in having slow muscle oxygen uptake dynamics close to that of pulmonary oxygen uptake. Several studies have considered whether the observed dynamics of oxygen uptake at the onset of exercise is the manifestation of an “inertia” in the rate of O2 delivery to the muscle fibers18 or of an intrinsic slowness of the intracellular oxidative metabolism.23,28

Linking Simulated and Measured Fluxes

The model developed here takes into account the dynamic interplay of oxygen delivery (via convection and diffusion),42 and mechanisms of cellular energy metabolism occurring in skeletal muscle during exercise.62 Specifically, muscle oxygen saturation (Fig. 2) and pulmonary oxygen uptake (Fig. 3) measurements were linked using a mathematical model to quantify the metabolic fluxes involved in cellular respiration in skeletal muscle during rhythmic voluntary contractions on a cycle ergometer.

According to the model of cellular energy metabolism proposed by Vicini and Kushmerick,62 φOxPhos is coupled to the metabolic fluxes of \({\phi_{\rm CK}^{\rm r}}\) and \({\phi_{\rm CK}^{\rm f}}\) (Fig. 4). This model describes how the dynamics of net PCr breakdown \({(\phi_{\rm CK}^{\rm r} -\phi_{\rm CK}^{\rm f})}\) results in a drop of C PCr to supply the amount of ATP (Fig. 5) that φOxPhos cannot synthesize rapidly enough (Fig. 4) to maintain ATP homeostasis at exercise onset (Fig.  5). Similarly, to Vicini and Kushmerick,62 we described cellular energy balance in skeletal muscle without including the glycolytic pathway as a source of ATP synthesis. Under certain conditions, such as maximal or supramaximal exercise, as well as at the onset of high intensity exercise, the glycolytic contribution to ATP synthesis may be significant. The results of our simulation at the onset of heavy and very heavy intensity exercise would most likely be different, were the glycolytic contribution be taken into account. Specifically, we speculate that the rate of phosphocreatine breakdown would be reduced in proportion to the net rate of lactate formation from glycolysis throughout the exercise period.

The simulated dynamic response of muscle oxygen utilization (indicated by φOxPhos) was always significantly faster than the response of pulmonary oxygen uptake at all intensities (Fig. 3).

The dynamic response \({\tau_{{\rm UO}_{\rm 2m}}}\) of simulated UO2m or φOxPhos was always faster than \({\tau_{{\rm VO}_{\rm 2p}}^{\rm j}}\) of measured VO2p at all exercise intensities \({(\tau_{{\rm VO}_{\rm 2p}}^{\rm j}}\), M: 65 ± 7 s; 100 ± 24 s; 82 ± 31 s).42 Furthermore, \({\tau_{{\rm UO}_{\rm 2m}}}\) based on model simulation in which φOxPhos depends on oxygen concentration Eq. (7) was 2–3 s larger than \({\tau_{{\rm UO}_{\rm 2m}}}\) obtained using an empirically derived exponential function for UO2m.42 This difference vanished when φOxPhos did not include oxygen dependence. The dependence of φOxPhos on oxygen concentration may become evident at low oxygen concentration values, especially at high-intensity exercise.

It should be noted that without the oxygen dependence, the mean V Vmax estimated is reduced from 44 to 37 mmol l−1 min−1. Therefore, the original expression for φOxPhos proposed by Vicini and Kushmerick62 was modified in the present model by including an oxygen dependency, which is important at higher exercise intensities where oxygen could be very low and affect metabolism.67–69 During exercise, oxygen dependency affects the dynamics of muscle oxygen utilization, but does not affect muscle oxygen uptake \({\tau_{{\rm VO}_{\rm 2m}}^{\rm j}}\) because the venous oxygen concentration dynamics is not greatly affected. It is worth noticing that, in our model, we assume the rate of ATP production associated with oxidative phosphorylation to be six times the rate of oxygen consumption (i.e., β = 6). Although the P/O2 ratio may be 5 or 6, as proposed in some authors,45,56 a recent review suggests values between 3 and 6.32 Nevertheless, when different values of β were used in the simulations, the estimated dynamics of oxidative phosphorylation remain the same at all exercise intensities, even though it affected the value of k jATPase and the estimation of V jmax .

The kinetic parameters of cellular metabolism Eq. (9–10) were previously estimated for nerve stimulation of finger and wrist flexors.62 The same parameter values were used in this study even though muscles engaged in bicycle exercise were voluntarily activated. Nevertheless, model simulations of the dynamics of oxidative phosphorylation and PCr, characterized by the mean response times \({\tau_{{\rm UO}_{\rm 2m}}}\) and τPCr, are in agreement with previously published studies performed in human subjects under similar conditions.53

The model equations account for cellular ATPase, oxidative phosphorylation and CK fluxes variations in skeletal muscle during exercise, where the cellular respiration is regulated by feedback control with dependence from ADP and oxygen concentrations.13–15,67–69 It is worth to note that, in the literature, several feedback and feed forward control models have been proposed to describe the regulation of cellular respiration in vivo, such as (1) feedback control using a Michaelis-Menten relationship between oxidative phosphorylation and [ADP], (2) higher-order feedback control using an expression in which the Hill coefficient is greater than 1, or (3) a more fundamental expression of relating oxidative phosphorylation to the free energy of ATP hydrolysis.34,41,63 Furthermore, some other scientists have proposed feed forward mechanisms to control oxidative phosphorylation.2,40 However, the experimental data available in the present study are not sufficient to address this issue; therefore, we adopted the approach successfully applied previously to in vivo data obtained by NMR spectroscopy.62

The computation of C ADP at steady-state during constant work rate exercise (M: 2.74 · 10−3 mM; H: 5.94 · 10−3 mM; V: 1.54 · 10−2 mM), based on the equilibrium assumption using the CK equilibrium constant K eq = 177 (pH = 7, T = 38 °C)57 and C ATP, C PCr, C Cr obtained from simulations using full kinetic expressions for CK (Eqs. 9–10), is consistent with the C ADP values simulated (M: 2.32 · 10−3 mM; H: 5.26 · 10−3 mM; V: 1.45 · 10−2 mM). Moreover, our approach allows for dynamic information on the metabolic fluxes, which is useful to investigate the mechanisms regulating the dynamics of oxidative phosphorylation during exercise.

Model Transport Parameter

In this model, mass transport of oxygen between capillary blood and tissue depends on permeability-surface area coefficient PS, which varies with muscle blood flow (Eq. 13) that in turn depends on exercise intensity.11 The dependency of PS on Q jm allows the model to simulate resting and exercise conditions. Although the expression used to describe the variation of PS with work rate is phenomenological, it is consistent with tissue oxygen concentration values at rest and in agreement with the hypothesis of a capillary recruitment during exercise.33 Regardless whether there is capillary recruitment during exercise or just an increase in blood flow rate through the already recruited capillaries, in our simulations, PS coefficient must be set to a sufficiently high value to ensure enough oxygen supply to contracting muscle fibers to match the energy demand at all intensities.42 The parameter Q C that affects the PS−Q jm relationship (Eq. 13), must be sufficiently small so that simulations of arterio-venous oxygen differences are consistent with the blood flow increase with exercise intensity. Also, from Eq. (12), the blood flow increases between steady states, ΔQ jm were in agreement with those found previously at all intensities.42

Limitations of Experimental Data

Model simulations indicate that a reduction of the amount of active muscle, V m, from 40 to 9 l and thus a reduction of capillary and tissue volume (V cap, V tis) engaged during exercise, affects the temporal profile of the metabolites, the maximum rate of oxidative phosphorylation V jmax and the mean response time of cellular oxygen utilization \({\tau_{{\rm UO}_{\rm 2m}}^{\rm j}}\). In particular, the smaller the V m values are the faster the dynamics of oxygen concentration in capillary blood and in tissue appears, for a fixed value of V jmax . The V jmax value estimated by Vicini and Kushmerick62 (30 mM min−1) falls in the range of values determined in our simulations (45–25 mM min−1) when varying V m between 40 and 9 l. The estimation of the V jmax parameter is limited by the experimental uncertainty in determining active muscle volume and the microvascular volume distribution during exercise. In the present study, the volume of muscle recruited by each subject during exercise was assumed to be a percentage of the whole body volume; thus, a more accurate determination of the functional relationship between active muscle volume engaged and work rate is needed. Experimental measurements of key metabolites, such as PCr in large muscle groups, could be coupled to the model developed to estimate the active muscle volume unknown. Future model development will take into consideration motor unit recruitment pattern accounting for an increase of muscle volume during exercise.

Another limitation related to the unknown muscle volume engaged during exercise is the lack of accuracy and specificity of NIRS measurements in contracting skeletal muscle. Specifically, the limitation of using NIRS technique to evaluate StO2m in skeletal muscle is its lack of discriminatory power to distinguish among the relative contributions of skin, adipose tissue,43 capillaries, and small arterioles and venules to the NIRS signal obtained from an uncertain volume of the region investigated. Therefore, the absolute values may be misleading with respect to the rate of adjustment of oxidative metabolism. Changes in capillary blood volume can affect the concentrations of oxygenated and deoxygenated hemoglobin measured by NIRS and cause variability in the evaluation of oxygen consumption.10,18 Another potential source of interference is the relative contribution of Hb and Mb44,60 to the NIRS signal. Most studies focus on Hb changes, since it has been reported that intracellular Mb accounts for less than 10% of the total NIRS signal.28 Currently, available NIRS instrumentation cannot accurately determine the relative contribution of myoglobin (Mb) to the total NIRS signal.

Future Directions

The complexity of the model required in future developments depends on the experimental information available and hypotheses to be tested. The validity of model simulations depends on the accuracy and precision associated with estimating unknown parameters. In modeling the linkage between whole body, tissue, and cellular processes with mass balances, the analysis of dynamic responses depends on the estimated volumes of the tissues.

In practice, if the observation scale is reduced to the microvascular level, the volume of perfused tissue under consideration is uncertain. Furthermore, the spatial distribution and temporal variation of blood flow and oxygen concentration can have a significant effect on the interpretation of the measurements.48 With appropriate experimental data, more general detailed models,5,6,17,22 can take into account (a) spatial distribution of transport and reaction processes in the capillaries and extra-vascular tissue, (b) intracellular compartmentation,46 and (c) structural and metabolic characteristics of the tissue (e.g., muscle fiber types). Specifically, for human studies of skeletal muscle metabolism in response to exercise, improvements in the analysis of pulmonary and muscle oxygen dynamics could be made by experiments that simultaneously combine invasive and non-invasive techniques.18,28,29,39,50,59 The variables to be measured are venous oxygen content, muscle blood flow, and muscle oxygenation during exercise. Measurements of muscle oxygen utilization should be as close as possible to the site of cellular respiration in order to assess the dynamics of oxidative phosphorylation in vivo during exercise. Measurements of 31P and 17O by magnetic resonance spectroscopy could provide key species concentrations in oxidative phosphorylation.16

Nowadays measurements of pulmonary oxygen uptake, VO2p are used to make inferences on the bio-energetic processes involved during exercise. The VO2p average measured in our population at the end of constant work rate exercise may be expressed as percent of the maximal pulmonary oxygen uptake.42: 49 ± 11%, 68 ± 9% and 87 ± 6% at moderate, heavy and very heavy exercise intensity, respectively.

Presence of the slow component was discernable from our data during heavy and very heavy intensity exercise. To date, there is no unequivocal explanation for such behavior because several factors can lead to an excess VO2p during high intensity constant work rate exercise, as already reported by several authors.4,20,21, 52 The predominant portion of the slow component (86%) may be associated with the contracting muscles,51 whereas fiber recruitment at higher level may affect the pattern of VO2p measured at the mouth. Fiber recruitment may only contribute in part to the slow component. The model developed here does not incorporate components that may lead to a slow component for VO2; however, it can be enhanced to enable testing of putative mechanisms resulting in excess oxygen consumption during high intensity exercise.

Additionally, if dynamic measurements of O2, CO2, and pH were obtained in the venous blood perfusing the contracting muscle, then the effect of acid-base regulation1,17 could be modeled to investigate the control of respiration during exercise. Although direct in vivo measurements of UO2m are not feasible during exercise involving a large muscle mass, indirect evaluation of UO2m can be achieved with a multi-scale systems biology approach. This would combine computational modeling and simulation with experimental data with many types of measurements at the whole body, tissue, and cellular levels.