## Abstract

On the basis of experimental studies, the intracellular O_{2} (iPo_{2})-work rate (WR) relationship in skeletal muscle is not unique. One study found that iPo_{2} reached a plateau at 60% of maximal WR, while another found that iPo_{2} decreased linearly at higher WR, inferring capillary permeability-surface area (*PS*) and blood-tissue O_{2} gradient, respectively, as alternative dominant factors for determining O_{2} diffusion changes during exercise. This relationship is affected by several factors, including O_{2} delivery and oxidative and glycolytic capacities of the muscle. In this study, these factors are examined using a mechanistic, mathematical model to analyze experimental data from contracting skeletal muscle and predict the effects of muscle contraction on O_{2} transport, glycogenolysis, and iPo_{2}. The model describes convection, O_{2} diffusion, and cellular metabolism, including anaerobic glycogenolysis. Consequently, the model simulates iPo_{2} in response to muscle contraction under a variety of experimental conditions. The model was validated by comparison of simulations of O_{2} uptake with corresponding experimental responses of electrically stimulated canine muscle under different O_{2} content, blood flow, and contraction intensities. The model allows hypothetical variation of *PS*, glycogenolytic capacity, and blood flow and predictions of the distinctive effects of these factors on the iPo_{2}-contraction intensity relationship in canine muscle. Although *PS* is the main factor regulating O_{2} diffusion rate, model simulations indicate that *PS* and O_{2} gradient have essential roles, depending on the specific conditions. Furthermore, the model predicts that different convection and diffusion patterns and metabolic factors may be responsible for different iPo_{2}-WR relationships in humans.

- oxygen diffusion
- convection
- bioenergetics
- exercise
- glycolysis
- contraction

the skeletal muscle intracellular Po_{2} (iPo_{2})-exercise intensity relationship was previously investigated to study the regulatory mechanisms of O_{2} consumption (V̇o_{2}) in working skeletal muscle (49, 55, 57). In particular, from these NMR studies in humans, skeletal muscle iPo_{2} was inferred from measurements of deoxymyoglobin concentration during muscle contraction. The relationships between iPo_{2} and metabolic intensity [i.e., work rates (WR)], however, were dissimilar. Molé et al. (49) reported that iPo_{2} decreased linearly as WR increased above 50–60% of peak V̇o_{2} (V̇o_{2peak}), whereas Richardson et al. (55, 57) found that iPo_{2} was low and constant above 60% of V̇o_{2peak} during knee extensor exercise. In addition, Molé et al. inferred that the blood-tissue O_{2} gradient, rather than capillary permeability-surface area (*PS*), as proposed by Richardson et al., was the dominant factor determining O_{2} diffusion. *PS* is an overall transport coefficient equal to the rate of O_{2} transport between capillary blood and myocytes divided by the corresponding O_{2} concentration difference. It is directly related to diffusion conductance or diffusing capacity. Previous mathematical models suggested that O_{2} diffusion primarily depends on the spacing and clearance of RBCs (18, 26), two factors that contribute to *PS*. Therefore, the role of the O_{2} gradient and *PS* on the regulation of O_{2} diffusion at different muscle contraction intensities is still not conclusive.

Because these studies were done on different muscles undergoing different exercise protocols, they are difficult to compare directly. V̇o_{2} and blood flow (Q̇) distribution in skeletal muscle likely differ in the two experimental protocols. In particular, single-legged knee extensor exercise results in much higher Q̇ and V̇o_{2} per unit muscle than does plantar flexion or two-legged cycle ergometry (6, 49, 56). Confounding the interpretation of these results further is the heterogeneity of Po_{2} within the muscle, which affects NMR measurements and subsequent estimates of iPo_{2}.

O_{2} delivery and cellular metabolism are the main processes that can limit the dynamic response of skeletal muscle V̇o_{2} during contraction (63). Therefore, besides O_{2} convection and diffusion effects on iPo_{2}, the mechanisms regulating the utilization of aerobic and anaerobic energy sources must be taken into account in quantifying the iPo_{2}-WR relationship during higher metabolic demand. Previous studies using the canine gastrocnemius model have provided insightful information on the factors limiting muscle oxidative metabolism during contraction (23–25). At submaximal metabolic intensity (23, 24) (e.g., 60% of V̇o_{2peak}), V̇o_{2} kinetics responses appear to be limited primarily by cellular metabolism; at higher metabolic intensity (100% of V̇o_{2peak}), O_{2} delivery (i.e., diffusion and convection) and metabolic processes appear to play a role (25). However, these studies of skeletal muscle contraction did not evaluate iPo_{2} or the rate of O_{2} diffusion.

In spontaneously perfused muscle, Q̇ increases with V̇o_{2} and energy demands of the muscle. Imposing an elevated Q̇ would not be expected to alter V̇o_{2} but might affect iPo_{2} or other key metabolites. Although the effects of Q̇ on the kinetics of V̇o_{2} have been determined experimentally, the effects of elevated Q̇ and *PS* on iPo_{2} in canine muscle are not known.

A quantitative approach proposed by Wagner (64) was successfully applied to understand the interaction between O_{2} convection and diffusion and their roles in determining maximal V̇o_{2}. However, this approach is limited to steady-state analysis, and the interaction between O_{2} delivery and aerobic and anaerobic systems is not considered. In response to muscle contraction, the metabolic fluxes associated with the glycolytic and oxidative systems are regulated by compensatory mechanisms involving ATP/ADP and NADH/NAD changes (7). Therefore, glycolytic and oxidative systems are indirectly related through the metabolic fluxes that share common substrates. Thus anaerobic glycolysis can indirectly affect iPo_{2} in skeletal muscle at higher exercise intensity.

To quantify the relationships between iPo_{2}, O_{2} utilization, and WR in skeletal muscle at higher exercise intensity, a mechanistic, mathematical model of O_{2} transport and metabolism can be used. For this purpose, a model of anaerobic glycogenolysis (44) was combined with a model of O_{2} transport and utilization (43). The resulting model delivers O_{2} by convection and diffusion and supplies energy by oxidative and glycogenolytic mechanisms. Accordingly, it can simulate skeletal muscle responses to submaximal and maximal metabolic demands. Model outputs compared with experimental data obtained under different experimental conditions provide model validation. Model simulations of canine skeletal muscle during contraction show how the iPo_{2}-energy demand relationship is affected by *PS*, blood perfusion patterns, and oxidative and glycogenolytic characteristics. Differences in O_{2} delivery and energy utilization patterns associated with various experimental conditions can explain and reconcile apparent differences in the iPo_{2}-WR relationship at different levels of metabolic demand. Since simultaneous measurements of tissue metabolites and iPo_{2} with contracting muscle are difficult to obtain, a mechanistic, mathematical model is used to simulate the metabolic response to higher metabolic rates. These simulations can provide a reliable basis for evaluating *1*) the importance of O_{2} gradient and *PS* in the control of O_{2} diffusion and *2*) the conditions under which iPo_{2} decreases linearly or remains constant in response to increasing muscle contraction intensity.

### Glossary

- C
_{A}^{T} - Total concentration of ATP, ADP, and AMP
- C
_{C}^{T} - Total concentration of phosphocreatine and creatine
- C
_{Hb} - Concentration of Hb in blood
- C
_{N}^{T} - Total concentration of NAD and NADH
- C
_{mc,Mb} - Concentration of myoglobin in myocytes
- C
_{rbc,Hb} - Concentration of Hb in RBCs
*D*_{O2}_{,b}- Effective dispersion coefficient of O
_{2}in blood *D*_{O2}_{,c}- Effective dispersion coefficient of O
_{2}in tissue - f
_{b} - Blood capillary volume fraction in muscle
- Hct
- Hematocrit
*k*_{ATPase}- ATPase rate constant
*K*_{i}- Rate coefficient for a reaction involving metabolite
*i* *K*_{Hb}- Hill constant at which Hb is 50% saturated by O
_{2} *K*_{Mb}- Hill constant at which myoglobin is 50% saturated by O
_{2} *M*_{lac}- Michaelis-Menten transport constant of lactate
*M*_{pyr}- Michaelis-Menten transport constant of pyruvate
*n*H- Hill coefficient
- OxP
- Oxidative phosphorylation
- Po
_{2,ven}^{R} - Resting Po
_{2}in the vein *PS*^{rest}- Permeability-surface area at rest
*PS*^{max}- Maximal value of permeability-surface area
- T
_{lac} - Maximal transport flux of lactate
- T
_{pyr} - Maximal transport flux of pyruvate
- V̇o
_{2}^{R} - Resting muscle V̇o
_{2}from experimental data - V̇o
_{2}^{S} - Steady-state exercise muscle V̇o
_{2}from experimental data *V*_{max}- Maximal flux rate
- V
_{mus} - Total muscle volume
- W
_{mc} - Volume fraction of myoglobin in myocytes
- α
_{gly} - Glycogenolytic capacity
- α
_{O2} - O
_{2}solubility in blood - γ
_{O2}_{,b} - Coefficient that accounts for the equilibrium between free and bound O
_{2}in blood - γ
_{O2}_{,c} - Coefficient that accounts for the equilibrium between free and bound O
_{2}in tissue

### Subscripts and Superscripts

- A-V
- Arteriovenous difference
- art
- Artery
- B
- Bound O
_{2}concentration - b
- Blood
- c
- Cell
- cap
- Capillary
- exp
- Experimental data
- F
- Free O
_{2}concentration - f
- Forward reaction
- mod
- Model simulation
- mus
- Muscle
- R
- Resting condition
- r
- Reverse reaction
- S
- Contraction condition
- SS
- Steady-state condition
- T
- Total O
_{2}concentration - tis
- Tissue
- ven
- Venous

### Enzymes

- ADK
- Adenylate kinase
- ALD
- Fructose bisphosphate aldolase
- ATPase
- Myosin ATPase
- CK
- Creatine kinase
- ENOL
- Enolase
- GAPDH
- Glyceraldehyde-3-phosphate dehydrogenase
- GP
- Glycogen phosphorylase
- LDH
- l-Lactate dehydrogenase
- PFK
- 6-Phosphofructokinase
- PGI
- Phosphoglucose isomerase
- PGK
- Phosphoglycerate kinase
- PGLM
- Phosphoglucomutase
- PGM
- Phosphoglycerate mutase
- PK
- Pyruvate kinase
- TPI
- Triose phosphate isomerase

### Metabolites

- ADP
- Adenosine diphosphate
- AMP
- Adenosine monophosphate
- ATP
- Adenosine triphosphate
- 13BPG
- Glycerate-1,3-biphosphate
- 2PG
- Glycerate-2-phosphate
- 3PG
- Glycerate-3-phosphate
- Cr
- Creatine
- DHAP
- 1,3-Dihydroxyacetone phosphate
- F6P
- Fructose-6-phosphate
- FBP
- Fructose-1,6-biphosphate
- G1P
- Glucose-1-phosphate
- G6P
- Glucose-6-phosphate
- GAP
- Glyceraldehyde-3-phosphate
- GLY
- Glycogen
- LAC
- l-Lactate
- NAD
- Nicotinamide adenine dinucleotide (oxidized)
- NADH
- Nicotinamide adenine dinucleotide (reduced)
- PCr
- Phosphocreatine
- PEP
- Phosphoenolpyruvate
- P
_{i} - Inorganic phosphate
- PYR
- Pyruvate

## METHODS

A mathematical model is developed to simulate and predict dynamic metabolic responses of skeletal muscle in vivo to a change in energy demand and experimental conditions. This model (Fig. 1) is based on two earlier models: *1*) O_{2} transport and metabolism in blood and tissue (43) and *2*) anaerobic glycogenolysis with key metabolites (including lactate and pyruvate) (44). The glycogenolytic pathway in the model specifically refers to glycogenolytic metabolism that results in lactate accumulation, i.e., an increase in intramuscular lactate concentration and lactate release into the blood. This model assumes that chemical species change continuously in space between the arterial input and venous output in blood and in tissue cells. The total muscle volume consists primarily of capillary blood and extravascular muscle tissue: V_{mus} = V_{cap} + V_{tis}.

### Spatially Distributed, Metabolite Dynamics

#### Concentration dynamics in blood.

Free O_{2} concentrations in the capillary blood [C_{O2,b}^{F}(ν,*t*)] and in muscle cells [C^{F}_{O2}_{,c}(ν,*t*)] depend on time (*t*) and tissue location, as indicated by the cumulative muscle volume (ν) from the arterial input (ν = 0) to the venous output (ν = V_{mus}). On the basis of a dynamic O_{2} balance in blood (43)
_{b} = V_{cap}/V_{mus}, *D*_{O2,b} is effective axial dispersion coefficient in blood, *J*_{O2}^{b,c} is transport flux of O_{2} between blood and tissue cells, and γ_{O2}_{,b} is a coefficient that accounts for the equilibrium between free and bound O_{2} in blood (see appendix b). A similar dynamic mass balance describes lactate and pyruvate concentrations

#### Concentration dynamics in extravascular tissue.

The dynamics of free O_{2} concentration in muscle cells of tissue can be described as
_{c} = V_{tis}/V_{mus} = 1 − f_{b}, *D*_{O2}_{,c} is the effective axial dispersion coefficient in blood, and γ_{O2}_{,c} is a coefficient that accounts for the equilibrium between free and bound O_{2} in muscle cells. For *j* species other than O_{2} (see appendix a), the dynamic mass balance in tissue has the following general form
*j* is
_{S↔P} is the reaction flux rate for which the species *j* is involved as substrate (S) or product (P) and β_{j,S↔P} is the corresponding stoichiometric coefficient.

#### Boundary and initial conditions.

For capillary blood, we assume that the input O_{2}, lactate, and pyruvate concentrations from arterial blood are known and that the output O_{2}, lactate, and pyruvate concentrations leaving the capillaries have negligible gradients in the direction of Q̇

#### Concentration constraints.

The concentrations of ATP, ADP, and AMP, phosphocreatine (PCr) and creatine (Cr), and NADH and NAD are related by mass conservation of adenosine (A), Cr (C), and NAD (N) within the cell. For constant total concentrations and volumes of the tissue

#### Blood-cell transport fluxes.

The capillary-tissue transport of O_{2} between blood and muscle cells depends on the *PS* and their free O_{2} concentrations
*PS* is assumed to change with the same dynamic pattern (43). For a step-up change of Q̇ with characteristic time constant τ_{Q̇}, we assume
*PS*^{max} is the maximal value and *PS*^{rest} is the value at resting, steady state, when *t* ≤ *t*_{0}. When Q̇ is pump-perfused (PP), *PS* is set at a constant *PS*^{max}.

For lactate and pyruvate (*j* = lac, pyr), the net transport fluxes *J*_{j}^{b,c} between blood and tissue depend on their concentrations (C_{j,b} and C_{j,c}) according to carrier-mediated transport (see Table 7) (14)

#### Metabolic fluxes.

The mathematical model describes the metabolic reactions within the cells (Fig. 1), including the glycogenolytic pathway, and ATPase, creatine kinase, and adenylate kinase reactions, as well as oxidative phosphorylation. Most of the reaction fluxes are described by ordered bi-bi reactions (see appendix b) (44)
_{a}, C_{b}, C_{p}, and C_{q} are substrates and products of the enzymatic reaction *i*, respectively, and *K*_{a}, *K*_{b}, *K*_{p}, and *K*_{q} are Michaelis-Menten constants. The maximal metabolic flux reverse and forward for each *i* reaction are calculated using the Haldane relation (44)
*V*_{max,f,i} values do not change with contraction intensity.

#### Energy supply and demand.

Fluxes related to energy supply and demand change depending on experimental conditions. The energy demand is associated with the ATPase reaction described by mass action (44)
_{gly}^{ATP} = ϕ_{PFK} + ϕ_{PGK} + ϕ_{PK}), oxidative phosphorylation (ϕ_{OxP}), creatine kinase (ϕ_{CK}), and adenylate kinase (ϕ_{ADK}). At steady state, flux of energy demand is balanced by the fluxes of energy supply
_{2} utilization. Several of these fluxes are represented distinctly. The creatine kinase reaction is a Cleland (random bi-bi) reaction (44)
_{i} according to the previous model developed (5) to simulate ATP production in mitochondria of cardiac muscle

#### V̇o_{2}.

For comparison with minimally invasive measurements, model outputs include V̇o_{2} and O_{2} saturation in blood and tissue. The rate of muscle V̇o_{2} per volume of muscle is
_{A-V}^{T}(*t*) is the arteriovenous concentration difference of total (free and bound) O_{2}. The equilibrium relationship between total and free O_{2}, e.g., arterial O_{2} (C_{O2,art}^{F}), is presented in appendix b (43).

### Simulation Strategy

This model is intended to simulate and predict the spatially distributed concentrations of O_{2}, lactate, pyruvate, Cr, PCr, AMP, and metabolites of the glycogenolytic pathway. Specifically, the model is applied to analyze underlying metabolic and transport processes based on responses of the canine gastrocnemius muscle in situ to tetanic stimuli (23–25) under a variety of experimental conditions: *1*) normoxia or hyperoxia, *2*) spontaneous self-perfused (SP) or controlled PP Q̇ patterns, and *3*) O_{2} utilization rate equivalent to 60% or 100% of peak, steady-state V̇o_{2} (V̇o_{2}^{ss}) (Table 1). Canine muscle was used for validation because of a readily available database containing direct measures of metabolic rate and many of its related components under different conditions: Q̇, muscle contraction frequency, and arterial O_{2} content (23–25).

In model simulations, the measured arterial concentration of free O_{2} (C_{O2, art}^{F}) corresponding to normoxia or hyperoxia was applied. The simulations also incorporated the measured Q̇(*t*). Simulations of the corresponding O_{2} utilization rate are based on equivalent estimation of *k*_{ATPase} (Table 2).

There are uncertainties about the glycogenolytic enzyme capacity used in this model. In particular, *V*_{max,f,i}, reported in the computational study (44) to simulate the whole glycogenolytic process in skeletal muscle, was adapted to describe only the anaerobic glycogenolysis in our model. To lower the anaerobic glycogenolytic capacity to physiological levels, a glycogenolytic capacity fraction 0 < α_{gly} < 1 was introduced. The parameter α_{gly} modifies the maximal flux rates (44) of key glycogenolytic reactions (Table 3). This fraction does not apply to the reactions of oxidative phosphorylation, creatine kinase, and adenylate kinase.

Model validation is based on the simulations of C^{T}_{A-V}(*t*) and V̇o_{2}(*t*) under all conditions with the same basic model parameter values. The model can also predict the dynamics of spatially averaged concentrations and metabolic fluxes in tissue. A general expression for any spatially averaged variable in this model with respect to tissue volume is as follows
_{2} (〈C_{ATP}〉 and 〈C^{F}_{O2}_{,c}〉) and fluxes of oxidative phosphorylation and ATPase (〈ϕ_{OxP}〉 and 〈ϕ_{ATPase}〉). Under steady-state conditions, the potential metabolic and transport limitations become apparent from several spatially averaged outputs. The average iPo_{2} 〈iPo_{2}〉, which is directly proportional to the average free O_{2} concentration in cells 〈C^{F}_{O2}_{,c}〉, reflects the O_{2} available in the tissue through transport from the blood. The O_{2} gradient, the difference between O_{2} concentration in capillaries and tissue 〈C^{F}_{O2}_{,b}〉 − 〈C^{F}_{O2}_{,c}〉, contributes to the O_{2} transport between blood and tissue. 〈C^{F}_{O2}_{,b}〉 and 〈C^{F}_{O2}_{,c}〉 are related to extracellular Po_{2} (ePo_{2}) and iPo_{2} by O_{2} solubility in blood and tissue, respectively. A constant average C_{ATP} 〈C_{ATP}〉 indicates that energy supply is meeting energy demand. The maximal flux fraction α_{gly} was increased to explore the effects of glycogenolytic capacity on 〈iPo_{2}〉 at higher ATP demand. Although there is no evidence that the full glycogenolytic capacity is used during the canine muscle contractions that were modeled, this condition may be relevant for intense human exercise. Also, to examine the effects of O_{2} transport limitations on the 〈iPo_{2}〉-WR relationship, physiologically normal and elevated values of blood perfusion and *PS* were simulated.

### Parameter Estimation and Numerical Solution

Values of most parameters in this model (Table 7) are available from previous studies: transport parameters of pyruvate and lactate (14), concentrations of arterial lactate and pyruvate (50), and parameters of the metabolic fluxes (Tables 3–6) (44, 63a). These parameter values apply to all experiments and are independent of the input conditions. The value of T_{lac} was chosen so that simulation of lactate output was the same as that measured by Grassi et al. (23) at submaximal metabolic intensity under constant Q̇. T_{pyr} was similarly determined from pyruvate uptake data (20).

In these simulations, most model parameter values are independent of the experimental conditions; therefore, their values are the same. Other parameter values for simulation depend on the experimental conditions (Tables 1 and 2). Values of muscle volume (V_{mus}) were obtained from experimental measurements of muscle weight and typical muscle density (23–25). Values of *k*_{ATPase} for resting and working muscle and for *V*_{max,OxP} were estimated by least-squares optimization (43) using measurements of arteriovenous O_{2} difference obtained under hyperoxic conditions. These experiments were conducted with constant PP Q̇ and injection of a drug (RSR-13) that causes a rightward shift in the O_{2}-Hb dissociation curve (24). In simulating the muscle response to contraction, values of *k*_{ATPase} for resting and working muscle were estimated to fit steady-state V̇o_{2} data for each experimental condition (23–25). The value of *V*_{max,OxP} was unchanged for simulations under all experimental conditions.

The initial model concentrations depend on the experimental conditions: arterial metabolite concentrations and Q̇ (Tables 1 and 2). These are obtained by numerical solution of the model equations using parameter values of Tables 3–7. The steady-state concentration values are the initial conditions at which the stimulus is applied.

For numerical solution of the model equations, the spatial derivatives are discretized to convert partial differential equations to ordinary differential equations (method of lines) using a well-developed code (DSS/2) (61). Consequently, the model consists of a set of ordinary differential-difference equations that constitute an initial-value problem solved using a robust algorithm for stiff systems (DLSODE) (27). The parameter estimation was obtained by minimizing the objective function represented by the difference between model simulation and experimental data. The objective function is minimized by numerical optimization using an adaptive, nonlinear algorithm (DN2FB; http://www.netlib.org) (16). Depending on the operating point in parameter space, changes of just 10% in the parameters of interest can produce significant differences in the simulated outputs.

To simulate the metabolic response to contraction under blood SP conditions, Q̇ and V̇o_{2} data of canine muscle under different contraction frequencies (22–25, 28, 31, 39, 66), were used (see Fig. 8*A*).

## RESULTS

### Parameter Estimation and Model Validation

By fitting simulated arteriovenous O_{2} differences to experimental data obtained under the conditions reported in Table 1 (hyperoxia + RSR-13, PP), optimal estimates were obtained for oxidative phosphorylation maximal flux rate (*V*_{max,OxP}; Table 3) as well as *k*_{ATPase} in resting (*k*_{ATPase,R}) and working (*k*_{ATPase,S}) muscle (Table 2). The model simulations produced dynamic outputs comparable to experimental data (Fig. 2) of Q̇, arteriovenous O_{2} concentration difference (C_{A-V}), and muscle uptake (V̇o_{2} = Q̇C_{A-V}). The validity of the mathematical model was tested by comparison of simulated muscle responses under various conditions (Tables 1 and 2) with corresponding experimental data (23–25). Values of *V*_{max,OxP} and other common parameters, estimated from previous studies, were assumed to be the same under all experimental conditions (Tables 3–7).

In response to energy demand corresponding to muscle contraction, the model simulations produced dynamic outputs comparable to experimental data (Figs. 3–5) for Q̇, C_{A-V}, and V̇o_{2}. Simulated responses of other variables are predictions for which corresponding experimental data are not available. Model simulation of V̇o_{2}, obtained during normoxia with spontaneous (SP) and PP Q̇, shows that the dynamic and steady-state responses are not noticeably affected by the muscle Q̇ pattern at 60% of V̇o_{2peak} (Fig. 3). Apparently, differences in Q̇ are compensated by changes in C_{A-V}. With constant perfusion and hyperoxia, C_{A-V} shows a delayed dynamic response with a larger slope (Fig. 4). The delayed dynamic response with a larger slope is still evident in the V̇o_{2} dynamics, but the overall V̇o_{2} kinetics remain unchanged. However, under more frequent muscle stimulation corresponding to 100% of V̇o_{2peak}, elevated muscle PP Q̇ results in a faster V̇o_{2} dynamic response (Fig. 5).

In addition to these outputs, the mathematical model predicts a myoglobin (Mb) saturation fraction of ∼0.95 at rest. Also, during the first 30 s of muscle contraction, the Mb saturation fraction and the iPo_{2} decrease rapidly, which contributes to an increase of the capillary-tissue O_{2} gradient.

### Factors Affecting Metabolic Response to Higher Energy Demand

#### Metabolic responses to contraction.

To simulate responses to higher rates of muscle contraction or ATP demand, equivalent higher values of *k*_{ATPase} were used. The model simulations reported in Figs. 6 and 7 were obtained using model parameters corresponding to 100% of V̇o_{2peak} with PP Q̇ of 101 ml·100 g^{−1}·min^{−1} through the transition from rest to contractions. The *k*_{ATPase} value predicted at 100% of V̇o_{2peak} (i.e., maximal *k*_{ATPase} in the canine gastrocnemius muscle preparation) is labeled *k*_{ATPase,max} (Figs. 6–8). As ATP demand (i.e., *k*_{ATPase}) increases during submaximal muscle contractions, the average O_{2} (〈ePo_{2}〉-〈iPo_{2}〉) gradient (Figs. 6*C* and 7*C*) and muscle V̇o_{2} (Figs. 6*B* and 7*B*) increase, while 〈iPo_{2}〉 decreases (Figs. 6*A* and 7*A*). 〈C_{ATP}〉 is constant over a wide range of submaximal energy demand but decreases at higher energy demand (Figs. 6*D* and 7*D*).

#### Effects of glycolytic capacity.

Changes in glycolytic capacity were simulated by changing the maximal metabolic flux rate parameters *V^*_{max,f,i} associated with the glycolytic pathway as determined by α_{gly}. Increased α_{gly} changes the steady-state relationships of O_{2} gradient, 〈iPo_{2}〉, and V̇o_{2} with contraction intensity (Fig. 6). For the same ATP demand (but before significant 〈C_{ATP}〉 decrease), an increase in α_{gly} lowers V̇o_{2} (Fig. 6*B*) and O_{2} gradient (Fig. 6*C*) but increases iPo_{2} (Fig. 6*A*).

#### Effects of O_{2} diffusion.

Changes in the O_{2} diffusion rate from capillaries to myocytes were simulated by changing the capillary *PS* via its maximum value (*PS*^{max}). For different *PS*^{max}, the steady-state relationships of iPo_{2}, V̇o_{2}, O_{2} gradient, and 〈C_{ATP}〉 vary with energy demand (Fig. 7). At low or high ATP demand, lower values of *PS*^{max} result in lower 〈iPo_{2}〉 and higher O_{2} gradient. At lower energy demand, V̇o_{2} and C_{ATP} are unaffected by *PS*^{max} (Fig. 7), while at higher energy demand, V̇o_{2} and C_{ATP} are reduced at lower *PS*^{max}. The energy demand at which C_{ATP} begins to decrease is lower with lower *PS*^{max} (Fig. 7).

#### Effects of convective transport.

Figures 6 and 7 show the effects of increased PP Q̇ at all contraction intensities. In addition, 〈iPo_{2}〉 was simulated for physiological values of SP Q̇ and V̇o_{2}. A variety of experimental data were used (22–25, 28, 31, 39, 66) to determine an empirical linear relationship between steady-state Q̇ and V̇o_{2} (Fig. 8*A*). The relationship was used in the model to estimate *k*_{ATPase} for several metabolic rates under physiological (spontaneous) perfused Q̇ to examine the effects of Q̇ on 〈iPo_{2}〉, 〈ePo_{2}〉-〈iPo_{2}〉 gradient, V̇o_{2}, and ATP. Optimal estimates for *k*_{ATPase} were obtained for steady-state Q̇ values, which are model inputs. Under spontaneous Q̇, V̇o_{2}, O_{2} gradient, and C_{ATP} are similar to values obtained with elevated (controlled) Q̇ (not shown). However, the simulated 〈iPo_{2}〉 is higher with elevated muscle perfusion than with normally perfused muscle at low contraction rates (Fig. 8*B*). At higher contraction rates, 〈iPo_{2}〉 is similar under both flow conditions.

At higher energy demand, as represented by a large *k*_{ATPase}, a V̇o_{2} plateau and decline in 〈C_{ATP}〉 appear (Figs. 6*B*, 6*D*, 7*B*, and 7*D*). This occurs when iPo_{2} decreases enough (Figs. 6*A* and 7*A*) to exert feedback control on the oxidative phosphorylation flux. In this situation, oxidative phosphorylation and anaerobic processes (e.g., glycogenolysis, creatine kinase, and adenylate kinase) cannot meet the energy demand without a reduction in 〈C_{ATP}〉. For the same muscle oxidative capacity (*V*_{max,OxP}), higher glycolytic capacity (i.e., α_{gly}) or *PS* (i.e., *PS*^{max}) allows the metabolic system to perform at higher energy demand without a decrease in 〈C_{ATP}〉 (Figs. 6 and 7).

#### Comparison with human experimental data.

According to the proportionality between *k*_{ATPase} and ATP demand, it is assumed that *k*_{ATPase} is proportional to WR [i.e., (*k*_{ATPase} − *k*_{ATPase,R})/(*k*_{ATPase,max} − *k*_{ATPase,R}) = WR/WR_{max}]. Therefore, the simulated iPo_{2} values in contracting canine muscle were compared with those obtained experimentally in humans (57, 49). At rest, simulated 〈iPo_{2}〉 was similar to the NMR-measured iPo_{2} in humans (58). Simulations were performed with a range of *PS*^{max} values, with WR increased until C_{ATP} decreased to 70% of its initial value. For each *PS*^{max}, WR_{max} in the experimental canine model corresponds to specific *k*_{ATPase,max} values at Q̇ = 100 ml·100 g^{−1}·min^{−1} according to the linear V̇o_{2}-Q̇ relationship (Fig. 8*A*). Under spontaneous Q̇, the effect of *PS*^{max} on 〈iPo_{2}〉 is similar to that observed with elevated blood perfusion (Fig. 9). Qualitatively, the model's prediction of the iPo_{2}-WR relationship is similar to that found by Molé et al. (49) and Richardson et al. (57). At low WR, the simulated iPo_{2} is higher than that observed during knee extensor exercise (57). At higher WR, the iPo_{2} values measured during knee extensor (57) and plantar flexion (49) exercises were close to those simulated at low *PS*^{max}.

## DISCUSSION

Experimental studies with human skeletal muscles show apparently contradictory results with respect to the change of iPo_{2} with increased energy demand (WR). One study (49) showed that iPo_{2} decreased linearly as WR increased, while others (55, 57) found that iPo_{2} decreased with increased WR to a plateau ∼60% of V̇o_{2peak}. Because these studies examined different muscles during different types of exercise, their results are difficult to compare directly. Furthermore, reliable estimates of iPo_{2} are very difficult to obtain because of heterogeneities within the investigated muscles. However, these different relationships can be explained by a quantitative analysis with our mechanistic, mathematical model that simulates O_{2} transport and cellular metabolism in canine skeletal muscle. The canine muscle preparation was chosen to validate the mathematical model and to quantitatively relate the model-simulated 〈iPo_{2}〉 changes to energy demand under conditions of normal (spontaneous) and elevated Q̇, as well as different O_{2} diffusion and metabolic characteristics, during contraction. In contrast to the experimental conditions in human studies, the canine preparation allows control of muscle Q̇ to obtain key information on factors affecting muscle V̇o_{2} kinetics. Model simulations show that the different experimental iPo_{2}-WR relationships are valid for different experimental conditions. The response of 〈iPo_{2}〉 to ATP demand depends on O_{2} convection, O_{2} diffusion, and metabolic capacity. In particular, the rate coefficient of *PS* appears to be the main controller of O_{2} diffusion. Model simulations suggest that iPo_{2} can decrease linearly for high *PS* [similar to results of Molé et al. (49)] or remain constant for low *PS* under higher muscle contraction intensity [similar to results of Richardson et al. (55, 57)] (Figs. 7 and 9).

### Model Characteristics, Simulations, and Predictions

The mathematical model used in this study combines a model of anaerobic glycogenolysis (44) with a model of O_{2} transport and utilization (43). To simulate the 〈iPo_{2}〉-WR relationship at higher ATP demand, the model incorporates O_{2} transport (i.e., convection and diffusion) with aerobic and anaerobic processes. The aerobic processes include oxidative phosphorylation, while the anaerobic processes include the creatine kinase reaction and glycogenolysis to sustain ATP supply to meet the demand. The glycogenolytic model (44) was modified to simulate experiments with perfused canine muscle in situ under different conditions using a scaling parameter α_{gly} for the maximal glycogenolytic fluxes. When α_{gly} = 0.2, the glycogenolytic fluxes involving ATP (Φ_{gly}^{ATP}) contribute ∼10% of the ATP required to meet moderate ATP demand. This corresponds to the percentage found for a canine muscle preparation during moderate contraction intensity (11, 12). However, if α_{gly} is set at >0.2, then the simulated tissue lactate concentration and glycolytic contributions to the energy demand are greater than those found under physiological conditions for canine muscle contraction (30).

Model simulations describe C_{A-V} and V̇o_{2} responses to different Q̇ and O_{2} input concentrations to energy demand (Figs. 2–5). These mimic experimental responses from electrically stimulated canine gastrocnemius muscle (23–25). While steady-state ATP demand must be estimated for each condition, the model is able to successfully predict V̇o_{2} and C_{A-V} dynamics without adjustment of any parameters. Simulations (Figs. 2–4) show that, in the canine gastrocnemius, convective transport does not limit oxidative metabolism at submaximal energy demand (23, 24) but could limit it at maximal energy demand (Fig. 5), e.g., 100% V̇o_{2peak} (25). At maximal contraction intensity, the rise time of V̇o_{2} dynamics in response to a step increase of energy demand is shorter under elevated PP (controlled) Q̇ than under SP (spontaneous) Q̇. This indicates a convective limit to O_{2} transport rate from blood to tissue, because controlled, elevated Q̇ can meet O_{2} demand more quickly. Nevertheless, factors such as oxidative and glycolytic capacity and O_{2} convection (SP and PP) patterns, as well as diffusion, can affect the relationship between iPo_{2} and energy demand. These results provide the basis for model validity with respect to mechanisms of the V̇o_{2} dynamics. Consequently, the model is sufficiently robust to predict underlying responses associated with O_{2} transport and metabolism that have not been measured.

### Relationship Between iPo_{2} and Energy Demand

#### Effect of PS on intracellular O_{2} and O_{2} diffusion.

When the O_{2} transport rate increases to match ATP demand, the O_{2} gradient adjusts according to *PS*^{max}. With increased energy demand from low to high contraction rate, the blood-tissue O_{2} gradient increases linearly (Fig. 7*C*), which causes O_{2} diffusion to increase linearly as well. At lower *PS*^{max}, the O_{2} gradient is increased for the same *k*_{ATPase} value (Fig. 7*C*) to provide enough muscle O_{2} delivery and sustain the same muscle V̇o_{2} and O_{2} utilization at submaximal energy demand (Fig. 7*B*). Consequently, at every metabolic rate, the O_{2} gradient is higher and the average iPo_{2} or 〈iPo_{2}〉 is lower than the values obtained at higher *PS*^{max} (Figs. 7*A* and 9), although V̇o_{2} is unchanged (Fig. 7*C*). Thus, under physiological conditions, *PS*^{max} can affect 〈ePo_{2}〉 and 〈iPo_{2}〉 temporal profiles during muscle contraction. At higher energy demand and lower *PS*^{max}, the O_{2} gradient and, therefore, O_{2} diffusion cannot increase further, because 〈iPo_{2}〉 is close to zero (Fig. 7*A*). As a result, oxidative phosphorylation is limited according to *Eq. 17* and V̇o_{2} reaches a plateau (Fig. 7*B*). Higher *PS*^{max} allows O_{2} diffusion and, therefore, V̇o_{2} to reach plateaus higher than those obtained at lower *PS*^{max}. At *PS*^{max} = 200 l·l^{−1}·min^{−1}, plateaus occur for V̇o_{2} = 16 ml O_{2}·100 g^{−1}·min^{−1} and for 〈iPo_{2}〉 = 10 mmHg (Fig. 7, *A* and *B*). Although 〈iPo_{2}〉 is greater than the critical Po_{2}, iPo_{2} is close to zero in some areas of the muscle tissue (simulations not reported), limiting O_{2} diffusion and utilization. Besides this transport limitation, there is also metabolic limitation resulting in a V̇o_{2} plateau.

These simulation results predict that while 〈iPo_{2}〉 can decrease at higher WR [consistent with some data (49)], *PS*^{max} is the major regulator of O_{2} diffusion between blood and tissue in skeletal muscle [consistent with other data (57)]. For the same energy demand, *PS*^{max} and the O_{2} gradient are not alternative factors regulating O_{2} diffusion. Rather, *PS*^{max} is the controlling factor that determines 〈iPo_{2}〉 and the O_{2} gradient. Even though *PS*^{max} regulates O_{2} diffusion, the O_{2} gradient is important under various physiological conditions. NMR studies (10, 49, 55, 57) report a rapid decline in Mb saturation fraction at the onset of contraction, which reflects a rapid decrease in iPo_{2}. This dynamic response is also evident in model simulations. As a result, the simulated O_{2} gradient from blood capillary to tissue increases, which enhances O_{2} diffusion during muscle contraction, as indicated by other modeling (43) and experimental (62) studies.

*PS* can also affect the relationship between 〈iPo_{2}〉 and WR (%WR_{max}; Fig. 9). Simulations show that 〈iPo_{2}〉 decreases significantly as %WR_{max} increases in the lower range of WR. The rate of decrease is lower at higher WR, especially for lower values of *PS*^{max}. These results indicate that the experimental iPo_{2}-%WR_{max} relationships reported in the plantar flexion (49) and knee extensor (57) studies, which appear to be different, may be produced by different muscle *PS*, which depends on the experimental conditions.

#### Resolution of apparent experimental discrepancies.

These simulations explain quantitatively how O_{2} diffusion may limit muscle oxidative metabolism at maximal contraction rate, as suggested from canine muscle data (25, 29). However, at submaximal contraction, V̇o_{2} on-kinetics are not affected by arterial O_{2} content, implying the absence of O_{2} diffusion limitation (24). Moreover, our model simulations indicate that the differences between the results of Molé et al. (49) and Richardson et al. (57) are not mutually exclusive. These differences can be explained by differences in the interaction of transport and metabolic processes between plantar flexion (49) and knee extensor exercise (57).

#### Effects of Q̇ on iPo_{2}.

In model simulations shown in Figs. 6 and 7, Q̇ is at a constant, elevated level of 101 ml·100 g^{−1}·min^{−1}. However, physiological Q̇ is often linearly related to V̇o_{2} (60). Therefore, model predictions of 〈iPo_{2}〉 under physiological Q̇ were investigated to see how they might differ from model predictions under elevated Q̇. At low contraction intensity, when spontaneous Q̇ is low, model predictions of the O_{2} gradient and 〈iPo_{2}〉 differ from those found when Q̇ is elevated (Fig. 8). At the same ATP demand, muscle utilizes the same amount of O_{2} to provide ATP, whether its Q̇ is spontaneous or elevated. At lower energy demand, O_{2} delivery is lower with spontaneous than with elevated Q̇. Contraction intensity has a greater effect on tissue O_{2} when Q̇ is low (Fig. 8*A*) than when it is elevated (Figs. 6 and 7). For the same V̇o_{2} and utilization at lower energy demand, steady-state 〈iPo_{2}〉 must be much lower in the spontaneous than in the elevated Q̇ condition. Simulations show that 〈iPo_{2}〉 at higher energy demand (Fig. 8*B*) is similar when spontaneous and elevated Q̇ are similar (Fig. 8*A*).

#### Effects of glycogenolytic capacity.

Greater glycolytic capacity (α_{gly}) is associated with higher 〈iPo_{2}〉, but lower V̇o_{2} (Fig. 6, *A–B*) and oxidative phosphorylation. Under this condition, glycogenolysis (through ϕ_{gly}^{ATP}) consumes more ADP and produces more ATP for the same energy demand. This lower ADP concentration reduces oxidative phosphorylation by feedback mechanisms; therefore, the oxidative contribution to the energy demand is smaller in muscle with higher glycolytic capacity. When α_{gly} is higher, 〈iPo_{2}〉 reaches a plateau at slightly higher energy demand due to the smaller oxidative contribution. In our mathematical model, 〈iPo_{2}〉 is indirectly related to the glycogenolytic capacity through the effect of ATP/ADP changes on the aerobic and anaerobic metabolic fluxes. These regulation mechanisms are in agreement with those proposed in a previous in silico study (7).

### Model Achievements, Limitations, and Future Developments

#### Factors affecting O_{2} transport.

Our model not only simulated an iPo_{2} plateau similar to the results of Richardson et al. (55), but it also predicted a V̇o_{2} plateau not yet observed. Model simulations show that, at the same energy demand and V̇o_{2}, iPo_{2} can be drastically different because of the capillary *PS* (Fig. 7). Corresponding to the experimental data (55), these simulations suggest that *PS*^{max} may increase with energy demand. In the mathematical model, *PS* is represented by an empirical time function associated with a step change in Q̇, which can be justified by experimental studies (8, 65). Although this model assumes that *PS*^{max} is constant and does not vary with energy demand, previous canine studies have reported an increase of *PS* with Q̇ and energy demand (19). Even when Q̇ is held constant, *PS* may increase with energy demand (17).

Changes in *PS* may be related to physical factors such as vasodilation, capillary recruitment, and RBC spacing in capillaries (18, 32, 52, 53), which have been proposed to explain an increase in O_{2} delivery during contraction. A mechanistic mathematical model is required to investigate the relationship between *PS* and these factors (18, 26) and their effects on O_{2} transport at the microvascular level. This study does not incorporate such mechanisms that determine *PS* changes but, rather, deals with the phenomenological role of *PS* and blood-tissue O_{2} gradient in determining the O_{2} flux. Regardless of the underlying mechanisms, *PS* increases at the onset of muscle contraction. The conclusions reached from this study are based on simulated *PS* changes during exercise, which are consistent with previous studies (18, 26, 55).

The mathematical model accounts for systemic Hct, but not capillary microvascular Hct. Capillary Hct, as measured in the microvasculature, is generally lower than systemic Hct due to differences in the velocity between plasma and RBCs (40). Our model was developed to provide a sufficiently detailed description of the overall O_{2} transport in the whole muscle without describing the microvascular effects. It allows direct comparison between simulated and measured Po_{2} in blood and tissue domains, as well as between simulated *PS* and O_{2} diffusion capacity experimentally determined by other groups (59). Other specific models and analysis are required to study the effect of Hct differences at the microvascular level (4).

The *PS* and O_{2} gradient of our model provide a phenomenological representation of O_{2} transport between RBCs and myocyte mitochondria. This model does not distinguish O_{2} gradients or *PS* within different blood and tissue domains. Consequently, model simulations are not able to distinguish, for example, cytosolic gradients of O_{2}. Such distinctions require a more detailed model (26). However, these distinctions are not expected to affect the overall relationships of *PS* and O_{2} gradient to iPo_{2}.

#### Factors affecting the iPo_{2}-WR relationship.

With our model, we were able to analyze factors affecting the iPo_{2}-ATP demand relationship in contracting canine muscle by simulating iPo_{2} in canine muscle at different WR. Comparisons are made to iPo_{2} measured in human muscle by NMR under knee extensor and plantar flexion (Fig. 9). This comparison, however, is limited by differences between humans and dogs, as well as differences in protocols, muscle fiber recruitment, maximal Q̇, and metabolic rate. Under the assumption that ATP demand changes proportionally with WR, the simulated results predict an iPo_{2}-WR relationship qualitatively similar to that observed in humans.

*PS* and the O_{2} gradient are not the only factors that can affect iPo_{2}. Other factors might include O_{2} transport between myocytes (34), Mb buffering within myocytes (33), and the critical Po_{2} of oxidative phosphorylation (45). The sensitivity of iPo_{2} to these factors may be relatively small, but important, under some conditions. Quantification of the significance of the model parameters *PS*^{max} and α_{gly} on iPo_{2} requires more experimental data under a wide range of physiological conditions.

Fiber recruitment may also affect the iPo_{2}-WR relationship in exercising human muscle (3, 47). The model applied in this work describes metabolic processes in electrically stimulated canine gastrocnemius muscle, where the predominant fibers are highly oxidative and all are recruited at the onset of muscle contraction. The simulated 〈iPo_{2}〉 reported in Fig. 9 is based on a mathematical model of canine muscle that does not take into account fiber-type composition. When type I and II fibers are recruited, muscle V̇o_{2} and metabolic control may differ from the case when only type I fibers are recruited. This is a consequence of differences in glycolytic and oxidative capacity and perhaps the P/O ratio of fiber types.

#### Interpretation of iPo_{2} changes with energy demand.

Our model simulations show that the iPo_{2} data of Richardson et al. (57) and Molé et al. (49) (Fig. 9) are not necessarily inconsistent when examined together. The extent to which differences exist between electrically stimulated canine gastrocnemius, human plantar flexion, and human knee extensor exercise may be attributed to experimental conditions that affect the maximal metabolic rate, which differs between muscles and species. The maximal power output is ∼6.6 times higher for knee extensor (54) than plantar flexion (49) exercise. This difference in power output seems to reflect differences in the metabolic state of each muscle at its WR_{max}; at 100% WR_{max}, ATP concentration remains at its resting value in plantar flexion (49) but is reduced by 30% in knee extension (2).

In isometrically stimulated canine gastrocnemius muscle, there is no power output; however, the maximal V̇o_{2} [16 ml·100 g^{−1}·min^{−1} (25)] is comparable to that in human plantar flexion [15 ml·100 g^{−1}·min^{−1} (49)] and lower than that in human knee extension [60 ml·100 g^{−1}·min^{−} (54)]. Maximal muscle perfusion also differs between the human knee extensor [385 ml·100 g^{−1}·min^{−1} (54)] and the isolated canine gastrocnemius in the present study [100 ml·100 g^{−1}·min^{−1} (25)]. However, the slope of the linear Q̇-V̇o_{2} relationship is similar for canine muscle (5.4 ml/ml; Fig. 8*A*), knee extensor [5.6 ml/ml (54)], and plantar flexion [6.1 ml/ml (6)]. Although there are protocol differences in maximal metabolic rate and Q̇, higher O_{2} delivery is balanced by higher V̇o_{2}. Therefore, these differences between experimental protocols are not expected to affect the iPo_{2}-WR relationship (Fig. 9). We note that V̇o_{2peak} and Q̇_{max} in canine perfused muscle (25) used to determine the maximal *k*_{ATPase} (Fig. 9; i.e., WR_{max}) are lower than those obtained for the same canine gastrocnemius model (V̇o_{2max} = 27 ml·100 g ^{−1}·min^{−1} and Q̇_{max} = 2.1 l/min) (38) with a different protocol. The V̇o_{2peak} commonly reported for canine gastrocnemius studies may not be the true V̇o_{2peak} achievable in this canine muscle preparation. However, addressing the limitations in reaching V̇o_{2peak} in a canine muscle preparation is beyond the scope of this report. If the maximal V̇o_{2} and Q̇ obtained by Kelley et al. (38) are used to determine %WR_{max} in the canine gastrocnemius, the model simulations reported in Fig. 9 would be leftward-shifted. However, these higher values do not alter the Q̇-V̇o_{2} relationship.

NMR-measured iPo_{2} (49, 55, 57) was compared with model-simulated 〈iPo_{2}〉 (Fig. 9). The correct interpretation of this comparison of simulated and measured iPo_{2} requires consideration of the NMR measurement limitations. In particular, iPo_{2} measurement is affected by P_{50} (i.e., Po_{2} at which 50% of Mb is bound to O_{2}). The P_{50} of skeletal muscle in vivo is still uncertain. Experimental values range from 1.5 to 5.5 mmHg (57), which can significantly change the quantitative measurement of iPo_{2} obtained by NMR studies. For comparison with the model simulations, the data of Richardson et al. (57) assumed P_{50} = 3.2 mmHg, while the data of Molé et al. (49) assumed P_{50} = 2.9 mmHg. For the model simulations, we assumed P_{50} = 2.42 mmHg. Therefore, quantitative comparison between simulated and measured iPo_{2} may be affected, because P_{50} may contribute to the different values of measured and simulated iPo_{2}. Although iPo_{2} calculation is affected by P_{50}, the iPo_{2}-WR relationship will have the hyperbolic trend reported in Fig. 9, regardless of the exact P_{50}.

#### Modeling oxidative and glycolytic systems.

The minimal mathematical model used in this study is sufficient for evaluating key factors related to O_{2} transport and metabolic limitations with respect to increased energy demand. The model components representing skeletal muscle glycogenolysis were adequate to analyze O_{2} diffusion limitation at higher energy demand. If sufficient experimental data, such as dynamic changes of glycogenolytic metabolite concentrations and/or fluxes with contraction intensity, were available, better estimates of model parameters could be obtained to further quantify mechanistic processes. Such data are needed to quantitatively analyze the potential mechanisms by which glycolysis is activated at higher WR (13, 21) and how glycogen limitations might affect V̇o_{2} (51). To account for pH-dependent enzyme kinetics, another model of glycogenolysis could also be incorporated (63a).

The simulated V̇o_{2} was obtained under the assumption that oxidative phosphorylation is controlled by ADP, P_{i}, and iPo_{2} according to previous experimental and modeling studies (5, 9, 43). In contrast, an experimental study (10) reported that V̇o_{2} is independent of ADP at the onset of contraction. For accurate prediction of key metabolite (e.g., ADP, NAD, and NADH) concentrations associated with regulation of phosphorylation and redox states in response to contraction, a model that distinguishes metabolic processes in cytosol and mitochondria is needed (46). A more detailed model of cellular respiration that includes oxidative phosphorylation, tricarboxylic acid cycle, and β-oxidation (5, 41, 46) would be worth incorporating if additional data were available to investigate regulatory mechanisms of energy metabolism in vivo (37). With this more detailed model, which would more accurately simulate NAD/NADH transport and metabolism, the relationship between iPo_{2} and the accumulation of anaerobic metabolic markers above the anaerobic threshold could also be investigated.

Experimental studies have observed a linear V̇o_{2}-WR relationship in skeletal muscle across a wide range of contraction intensities (1, 54, 60). Computational studies have attributed this linear relationship at higher exercise intensity to a parallel activation mechanism affecting oxidative phosphorylation (48). Our model is based on feedback control of ADP and P_{i} without the inclusion of detailed activation mechanisms previously proposed (48). Even so, it is able to successfully simulate V̇o_{2} kinetics in the canine gastrocnemius at 60% and 100% of V̇o_{2peak}. At submaximal (60% of V̇o_{2peak}; i.e., below *k*_{ATPase,max}) and maximal (100% of V̇o_{2peak}) contraction intensity (i.e., at *k*_{ATPase,max}), the model predicts a linear V̇o_{2}-ATP demand relationship (Figs. 6*B* and 7*B*). Above *k*_{ATPase,max}, this relationship reaches a plateau where V̇o_{2} is no longer linearly related to ATP demand (Figs. 6*B* and 7*B*). Human studies (15, 35) reported a plateau in the V̇o_{2}-WR relationship near maximal exercise intensity. However, uncertainties remain with regard to the maximal V̇o_{2} attainable in the canine muscle preparation (38). Because maximal V̇o_{2} at lower *PS*^{max} is primarily limited by O_{2} transport, inclusion of a parallel activation mechanism in our model would not significantly alter the simulated results and conclusion. Furthermore, in electrically (as opposed to neurally) stimulated muscle at the contraction intensities considered in this study, parallel activation may not occur (36, 42).

### Perspectives and Significance

Model simulations demonstrate that, depending on experimental conditions, the interaction of O_{2} diffusion with oxidative and glycolytic processes can lead to different relationships between iPo_{2}, V̇o_{2}, and ATP demand (49, 57) in the canine gastrocnemius. The results of model simulations indicate that *1*) *PS* is the main factor that regulates O_{2} diffusion and determines the O_{2} gradient during muscle contraction and *2*) iPo_{2} may decrease linearly or remain constant at high-intensity exercise according to different muscle O_{2} transport and metabolic responses to contraction. Consequently, apparently contradictory experimental results (49, 57) relating iPo_{2} and WR or metabolic rate in humans can be quantitatively reconciled.

Model simulations can predict quantitatively how O_{2} transport from capillaries to myocytes affects V̇o_{2} under different physiological and pathophysiological conditions. These possible variations can change not only steady-state V̇o_{2}, but also V̇o_{2} kinetics. In conjunction with experimental studies, model simulations provide a framework in which to distinguish transport and metabolic limitations on V̇o_{2} and iPo_{2}. Furthermore, model predictions provide an efficient basis on which to design critical, new experiments to investigate the distinctive effects of O_{2} convection and diffusion in limiting O_{2} utilization during muscle contraction.

## GRANTS

This research was supported in part by National Institutes of Health Grants GM-088823 and K25 AR-057206 and National Science Foundation Grant 0743705. J. Spires was supported by National Institute of General Medical Sciences Predoctoral Fellowship F31 GM-084682.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

J.S., G.M.S., and N.L. are responsible for conception and design of the research; J.S., L.B.G., B.G., G.M.S., and N.L. analyzed the data; J.S., L.B.G., B.G., G.M.S., and N.L. interpreted the results of the experiments; J.S. prepared the figures; J.S., G.M.S., and N.L. drafted the manuscript; J.S., L.B.G., B.G., G.M.S., and N.L. edited and revised the manuscript; J.S., L.B.G., B.G., G.M.S., and N.L. approved the final version of the manuscript; L.B.G. and B.G. performed the experiments.

## APPENDIX A

### Dynamic Mass Balance

The model is spatially distributed in blood and tissue compartments. All metabolites are present in the tissue compartment. The blood compartment only includes O_{2}, lactate, and pyruvate.

#### Blood compartment.

#### Tissue compartment.

### APPENDIX B

#### O_{2} Equilibrium

O_{2} can be free or bound to Hb or Mb. The sum of free and bound O_{2} is total O_{2} concentration
_{2} concentrations are related by local chemical equilibria. In blood, the relationship is

#### Transport Fluxes

The axial fluxes for the substrates (O_{2}, lactate, and pyruvate) present in the blood and tissue compartments are as follows

#### Metabolic Fluxes

Most of these fluxes are calculated using ordered bi-bi Michaelis-Menten kinetics. The exceptions are the creatine kinase flux, the oxidative phosphorylation flux, and the ATPase flux.

The metabolic flux expression for each reaction is reported as follows.

Glycogen phosphorylase
_{GP}) accounts for isozyme form *A* and *B* with ω_{a} and ω_{b} fractions as follows
*a* and *b*) to the other during contraction.

Phosphoglucomutase

Phosphoglucoisomerase

Phosphofructokinase
*L* are computed with the following expressions

Aldolase

Triose phosphate isomerase

Glyceraldehyde-3-phosphate dehydrogenase
*D* is calculated as

Phosphoglycerate kinase

Phosphoglyceromutase

Enolase

Pyruvate kinase

Lactate dehydrogenase

Adenylate kinase

- Copyright © 2012 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 63a.↵
- 64.↵
- 65.↵
- 66.↵