## Abstract

Cellular energy balance requires that the physiological demands by ATP-utilizing functions be matched by ATP synthesis to sustain muscle activity. We devised a new method of analysis of these processes in data from single individuals. Our approach is based on the logic of current information on the major mechanisms involved in this energy balance and can quantify not directly measurable parameters that govern those mechanisms. We use a mathematical model that simulates by ordinary, nonlinear differential equations three components of cellular bioenergetics (cellular ATP flux, mitochondrial oxidative phosphorylation, and creatine kinase kinetics). We incorporate data under resting conditions, during the transition toward a steady state of stimulation and during the transition during recovery back to the original resting state. Making use of prior information about the kinetic parameters, we fitted the model to previously published dynamic phosphocreatine (PCr) and inorganic phosphate (P_{i}) data obtained in normal subjects with an activity-recovery protocol using^{31}P nuclear magnetic resonance spectroscopy. The experiment consisted of a baseline phase, an ischemic phase (during which muscle stimulation and PCr utilization occurred), and an aerobic recovery phase. The model described satisfactorily the kinetics of the changes in PCr and P_{i} and allowed estimation of the maximal velocity of oxidative phosphorylation and of the net ATP flux in individuals both at rest and during stimulation. This work lays the foundation for a quantitative, model-based approach to the study of in vivo muscle energy balance in intact muscle systems, including human muscle.

- adenosine triphosphatase
- ischemia
- creatine kinase
- oxidative phosphorylation
- skeletal muscle
- SAAM II
- parameter estimation

the production of cellular work (electrical, chemical, or mechanical) must be coupled to energetically favorable splitting of ATP (ATPase). Oxidative phosphorylation (OxPhos) supplies the majority of the required energy. In muscle activity, this supply of chemical energy varies greatly from rest to exercise, providing data for investigating intracellular regulation of cellular respiration. Meyer (31) developed a model for this supply-demand matching in a feedback system based on an analogy with an electrical circuit. That model aided the physiological interpretation of the chemical changes, since it described experimental data well with a parsimonious explanation. Several other approaches exist, each with different features, to model (sometimes mechanistically) energy balance (10,15)

Several important elements are missing from these approaches. These models do not integrate kinetic expressions for relevant physiological mechanisms (oxidative phosphorylation, ATPase, and creatine kinase) and thus cannot estimate the kinetic and energetic parameters of those mechanisms from data. It is usually assumed instead that in vitro properties apply to the intracellular environment in vivo, because many properties are inaccessible to direct measurement in vivo and only composite quantities, such as the rate of change in phosphocreatine (PCr) concentration, are directly measured. Second, models usually only consider transitions from resting to activated state or the recovery phase of an exercise; the exception is Meyer's model (31), which contains a symmetrical constraint on the time courses during the onset and recovery from exercise. Our work will consider simultaneous modeling of both the resting state and the activated state, without constraints on the time course of exercise and recovery. A mechanistic and robust model coupled with data analysis could thus provide the needed noninvasive parameter estimation results for specific mechanisms included in the model. In human studies, such model analysis would be useful to repeatedly evaluate, in the same individual, responses to interventions such as adaptation to training or therapy. ^{31}P nuclear magnetic resonance (NMR) spectroscopy has opened new frontiers for the noninvasive quantification of cellular energy balance, making this kind of investigation possible in humans and providing a rich set of relevant measurements with good resolution (for reviews see Refs. 9, 33, 36). The quantities measured by NMR are the intracellular concentrations of high-energy phosphates, ATP and PCr, and of pH, all of which directly interact via the creatine kinase reaction
Equation 1The dynamics of PCr content in muscle cells (e.g., a decrease with increasing mechanical activity) provide a measure of intracellular ATP breakdown and recovery rates, with the assumption that the reaction in*Eq. 1
* remains close to equilibrium. The last assumption is not used in the present work; other analyses (25,33) show that equilibration is not critical for analysis. The muscle content of P_{i}, which also changes with activity (usually negatively proportional to PCr), is also quantifiable via NMR.

We define in this work an approach to the quantification of mechanisms essential to account for energy balance from NMR measurements made in the skeletal muscle of healthy human volunteers. The mathematical model consists of a minimal number of components, to assess the extent to which these well-defined and well-understood cellular and biochemical mechanisms may account quantitatively for the transitions from rest to exercise and back to rest. Thus our model contains elements and feedback concepts used previously (10, 31). If successful, such an analysis would provide an efficient method to estimate physiological and biochemical quantities that appear as parameters in the mechanistic equations of the model. These quantities help characterize energetics of muscle performance at a biochemical level and quantify changes in normal and pathological adaptations in human health and disease. Discrepancies between experimental observations and the model would indicate fruitful research directions.

We develop and illustrate our analysis with published data obtained from human muscle (6). Our goal was to match our model, with prior information on some of its parameters, to experimental data describing steady-state and transient changes in PCr and P_{i}. An important novelty of our work is that we analyze the entire time course of the experiment, that is, resting state to activity and recovery back to the resting state. We expected that the model would quantify some basic features of muscle cellular energetics, would confirm estimates made by empirical approaches (6), and would possibly allow additional parameters to be estimated. An important feature of the parameter estimation techniques employed is that the precision of the estimates correctly incorporates experimental error, enabling evaluation of their reliability in each experiment. Note, however, that such reliability is determined by the precision of the estimates and not by their accuracy, which would require independent model validation as discussed later. This report demonstrates our initial progress toward making this analysis. Once we proceed to the validation stage of model development, we expect to be able to use the model and subsequent embellishments of it to quantify in single individuals the fundamental biochemical characteristics of energy balance as specified by the equations for ATPase and ATP synthesis.

## EXPERIMENTAL DESIGN AND METHODS

### Data and Protocol

In this study, we used data already published in Ref. 6, to which we refer the interested reader for details. Table1 and Fig.1 summarize the phases of this experimental protocol and the expected conditions of ATPase flux and OxPhos flux. Cellular metabolite concentrations during nerve stimulation and muscle activation and subsequent recovery were measured in eight normal subjects using ^{31}P NMR.

#### Phase A.

From the start of the experiment to 180 s afterward (*phase A*, control period), control data were obtained during rest with normal blood flow to establish the baseline conditions for NMR data acquisition. The ATPase flux is at its basal rate in inactive muscle and is matched by OxPhos flux. All timing is nominal (seediscussion).

#### Phase B.

The delivery of oxygen to the muscle was stopped with a pressure cuff at 180 s (beginning *phase B*, resting ischemia). The muscle consumed stored oxygen at a low but finite resting rate, requiring ∼300 s of basal metabolism to entirely deplete the reserve. The muscle became completely anoxic by ∼480 s, staying so until 540 s. Verification of anoxia was given in the original report (6); after ∼5 min of ischemia, PCr was found to decrease while the muscle remained at rest, and the rate of decrease agreed with independent measures of resting metabolism. OxPhos flux stopped at anoxia, but ATPase flux continued at basal steady state.

#### Phase C.

From 540 to 668 s (*phase C*, anoxic muscle stimulation and PCr utilization), nerve activation with twitch stimulations every second stimulated the muscle, decreasing PCr and increasing P_{i}. Anoxia prevented OxPhos flux, but glycolytic ATP synthesis was possible.

#### Phase D.

From 668 to 757 s (*phase D*), stimulation stopped and the muscle relaxed while remaining anoxic. ATPase flux returned to basal; the curve as drawn in Fig. 1 assumes that glycolytic resynthesis matched the ATPase flux.

#### Phase E.

At 757 s, removal of the pressure cuff allowed blood flow to return to the arm from 757 to ∼1,400 s (*phase E*, aerobic recovery). OxPhos flux resynthesized PCr, restoring it to the initial concentration of resting muscle.

#### Phase F.

The original data (6) showed that the metabolite concentrations from ∼1,400 s onward (*phase F*, final resting state) were equal to those in *phase A*.

The original report (6) calculated absolute concentrations of P_{i} (mM) and PCr (mM) from direct NMR measurements using metabolite levels reported for muscle biopsy (18). These experiments also measured the time course of pH. All metabolite concentrations are in millimoles per liter of cell water, using the factor 0.7 as the fraction of intracellular water per kilogram of fat-free muscle (40).

### Mathematical Model

Kushmerick (25) has already partially described the model; it was used to simulate cell energy balance with approximate values for the relevant quantities. The equations used were based on established concepts of feedback control (e.g., Ref. 31) but full and realistic kinetic expressions were used. The unique feature of the model used here is that a common set of rate equations is used and constrained by the entire data set, including the resting, active, and recovering states.

ADP is the central regulated metabolite in this model, controlled by the balance between ATPase and ATP synthesis. Let us define the function ΔADP(*t*) as the difference between the sum of the initial (steady-state) ADP (ADP_{0}) and ATP (ATP_{0}) concentrations and the ATP concentration at*time t*
Equation 2Three components affect changes in ATP: *1*) decrease by ATPase flux, *2*) increase by OxPhos flux, and*3*) increase or decrease by the CK flux depending only on its substrate and product concentrations.

The following equation describes the change in ATP concentration over time
Equation 3where PCr_{0} is the initial value of PCr; note that the resulting fluxes will all be expressed as millimoles of ATP per liter, per second. The magnitude of the term dATP(*t*)/d*t* is small, because of the last term of*Eq. 3
*, which represents the temporal and spatial buffering of the ATP/ADP ratio in the cellular milieu (*Eq. 5
*). The parameter *k* (expressed in units of inverse time) is the rate constant for intracellular ATPase flux. It varies as a step function from its resting state value, *k*
_{rest}, to a different stimulated value, *k*
_{stim}, at the onset of muscle stimulation, returning to its initial value at the end of stimulation (or to other values, defined later). Because of the structure of the model, it is apparent that the parameter *k*represents the net balance of the specified processes plus any other ATP-utilizing and ATP-synthesizing pathway not explicitly stated in the differential equations.

A Hill sigmoid function describes the time course of the flux of ATP synthesis by OxPhos flux (21)
Equation 4where *V*
_{max} is the maximal velocity of OxPhos flux at infinite ADP concentration, *K*
_{ADP}is the substrate concentration at which the flux of ATP synthesis by OxPhos flux equals half of *V*
_{max} and*n*H is the Hill coefficient.

The equation describing PCr change over time is
Equation 5The constants *K*
_{ia},*K*
_{iq}, *K*
_{ib},*K*
_{b}, and *K*
_{p} are creatine kinase kinetic constants (30, 38). We did not account for the so-called “dead-end complexes” of the enzyme, because preliminary analyses showed the inclusion of appropriate terms did not alter the results and only made computation slower.*V*
_{for CK} is the maximal velocity in the forward CK flux, and *V*
_{rev CK} is the maximal velocity in the reverse CK flux. In the numerator, the first term indicates the flux in the direction of ATP formation, and the second term indicates the direction of PCr synthesis. The terms in the denominator denote the various forms of substrate and product binding to the enzyme.

The following equation describes the time course of P_{i}during PCr and ATP transient
Equation 6where P_{i0} is the steady-state value at rest for P_{i}. We will refer to the model defined by *Eqs.3
*–*
6
* as “*Model 1*” in the rest of this report.

### Parameter Estimation

Predicted PCr concentration, PCr(*t*), the solution to*Eq. 5
*, and P_{i} concentration, P_{i}(*t*), the solution to *Eq. 6
*, are fitted to data. To estimate the parameters of a generic mathematical model from data, a function of the difference between model prediction and data (the objective function) is minimized with respect to the model parameters (8). The classic approach to parameter estimation calls for the minimization of the weighted residuals sum of squares (WRSS)
Equation 7where *N* is the number of data points,*y _{i}
* is the data point at

*time t*, p̂ is the estimated parameter vector (size

_{i}*Np*), m(p̂,

*t*) is the model function evaluated at the parameter vector p̂ and

_{i}*t*, and V

_{i}_{i}is the measurement error variance of the data point

*y*. Minimization of this function allows the determination both of the optimal parameter values for each set of experimental data and of their precision. A lower bound on the covariance matrix (Cov) of the parameter estimates can be calculated from the inverse of the Fisher information matrix (8), which for uncorrelated errors is Equation 8where

_{i}^{T}indicates transpose. The diagonal of the covariance gives the variance of the estimates. The precision is usually expressed as percent coefficient of variation (%CV), i.e., the square root of the variances divided by the estimates' values.

The assumption behind *Eq. 7
* is that all model parameters are identifiable a priori from data; that is, all the information about the parameters comes from the kinetic data (12). However, as we have seen in our case, the model parameterization is too rich to ensure a priori identifiability of all parameters of interest. Fortunately, some of the parameters have independent information. Typically, this a priori knowledge consists of means and standard deviations obtained from population estimates. To fully incorporate this information, the objective function needs augmentation with an appropriate term (often called the “Bayesian” term), which tends to maintain the estimated parameter “close” to its independently measured value depending on the degree of constraint given by the standard deviation of that estimate
Equation 9where *N*
_{b} is the number of (uncorrelated) Bayesian parameters (smaller than or equal to the total number of parameters in the model *N*
_{p}), μ_{i} is the mean, Σ_{i} is the variance associated with the independent measure of the parameter*p _{i}
*, and p̂ is the estimated value of the parameter

*p*. This augmented objective function is called the maximum a posteriori (MAP) estimator objective function. A recent application in physiological modeling is in Ref. 11.

_{i}After defining the MAP objective function, we performed parameter estimation on the entire set of parameters *N*
_{p}. Note that it is still possible to calculate the precision of the parameter estimates using an extension of *Eq. 8
*
Equation 10With this augmented objective function, the fitting algorithm determines the optimal value of the parameters using both the data and any available prior information.

Experimental evidence suggested that the measurement errors affecting the P_{i} and PCr data were assumed to be Gaussian, zero mean, and with a constant standard deviation (P_{i}) and a constant coefficient of variation (PCr). Because of the method of obtaining the NMR data, these assumptions are not strictly correct throughout the time course of the experiment. The signal-to-noise ratio of the PCr peak integrals declines when PCr declines because of a fixed noise base; similarly, the signal-to-noise ratio of the P_{i} peak integrals increases as P_{i} increases. Thus the measurement error will change during the different phases of the experiment as PCr and P_{i} change. We have not yet accounted for such features of the data in our modeling. However, the magnitudes of the measurement errors were estimated a posteriori from the kinetic curves using extended least squares, as described in Ref. 5. For concision, we do not report here the exact definition of the objective function, which results in a further extension of *Eq. 9
*, or details on the procedure to estimate data error (4). The software SAAM II, version 1.1.1 (SAAM Institute and University of Washington, Seattle, WA) accomplished the model fitting, following the description of its kinetic analysis of metabolic and transient data in (4). SAAM II used the Rosenbrock differential equation integrator, with a relative accuracy of 0.1%, and established convergence of parameter estimation when objective function values did not change by more than 0.01% between iterations.

### Reduction of the Number of Adjustable Model Parameters

Some of the model parameters listed in Table2 are well known from a priori information. ATP_{0} and total creatine (TCr) have been measured in human muscle fibers and are independent of the cell type (18). Thus the initial condition for ATP is 8.2 mM, and TCr is 42 mM. The initial condition for PCr(*t*) (*Eq.5
*), PCr_{0}, and the value for ADP_{0} can both be uniquely calculated as a function of other model parameters from the steady-state equations
Equation 11and
Equation 12Estimates of the forward creatine kinase activity vary. We used 30 mM/s, a value less than the range of values reported for human muscle (100–200 mM/s; see Refs. 41, 43), because we did not use terms for the “dead-end” complexes, inclusion of which significantly reduces the effective enzyme activity by about fivefold (30). We established (not shown) that the analyses are not sensitive to this value in the range of 10 to 100 mM/s. The Haldane relationship constrained the maximal reverse fluxes as done previously (25)
Equation 13The initial pH is typically 7.0 to 7.1 in human muscle. The following relationship relates *K*
_{eq} of the creatine kinase reaction to pH measurements obtained during the course of the experiment
Equation 14

### Use of Available Prior Information as Bayesian Variables

Table 2 lists all the model parameters, several of which have been measured for human forearm muscle; they were fit as adjustable Bayesian parameters (see *Parameter Estimation*, above):*K*
_{ADP}, *n*H, and*k*
_{rest}. Prior mean and standard deviation values were taken from Ref. 21 for *K*
_{ADP} = 0.043 ± 0.003 mM and for *n*H = 2.11 ± 0.34. Mean and standard deviation for *k*
_{rest}(6) were 0.008 ± 0.001 per second for [ATP] = 8.2 mM; our present model also estimated basal ATPase rate as described below in *Independent Measure of Basal Rates of ATP Utilization. V*
_{max} and*k*
_{stim} were unconstrained, adjustable parameters with no prior information.

### Independent Measure of Basal Rates of ATP Utilization

Other experimental data on four of the experimental subjects (see Fig. 3 from Ref. 6) were available from our archive for analysis with our present model. The experiment consisted of an ischemic phase (corresponding to *A* in Fig. 1) lasting 1,200 s. After ∼500 s, PCr declined at the basal ATPase flux for ∼700 s more, corresponding to an extension of *phase B* in Fig. 1. PCr concentration declined measurably in the absence of stimulation.

Figure 2 shows a typical fit; the other three subjects gave equally good fits. Table3 displays the numerical results of this analysis. These values for *k*
_{rest} yield an average basal ATPase of 0.007 ± 0.001 mM/s (0.00087 s^{−1} · 8.2 mM), consistent with the result given in Ref. 6, cited above. This analysis of basal ATPase flux can constrain*k*
_{rest} to follow a certain mean and standard deviation. Therefore, we used a Bayesian prior on*k*
_{rest} equal to 0.00087 ± 0.00014 s^{−1}.

### Statistical Methods

Results are expressed as means ± SD. Differences between parameter estimates were evaluated using unpaired *t*-tests with a significance level set at 0.05. Model comparison was made through the Akaike information criterion (AIC; Ref. 1) and the Schwarz-Bayesian information criterion (BIC; Ref. 39). Briefly, these criteria weigh together complexity and goodness of fit of the models, by combining the regression objective function and the number of model parameters; between two models, the most parsimonious is the one that yields the lower criterion value. Akaike originally used AIC as a criterion for linear models; it has been claimed (3,37) that AIC applies to almost all statistical problems.

## RESULTS

### Model 1: Basic Model Analysis of Energy Balance Parameters

We begin our analysis with the conditions stated in Table 2; these are similar to the conditions used in prior simulations (25). Table 4 presents the results from the analysis of the eight subjects discussed here. The initial PCr concentration, PCr_{0}, and the initial ADP concentration, ADP_{0}, were all in the expected physiological range. On average, these were 0.016 ± 0.008 mM for ADP_{0} and 32.1 ± 3.7 mM for PCr_{0}. The estimated inorganic phosphate initial concentration, P_{i0}, averaged 3.17 ± 0.71 mM.

Table 4 contains all of the estimated parameters, and Fig.3 shows both the average of the model fits for the individuals and the average of the individual data. The model provided good estimates for the parameters for OxPhos flux. The average value for *V*
_{max} from these transient experiments is 0.50 ± 0.34 mM ATP/s, somewhat larger than the estimates for *V*
_{max} made from the maximal sustainable steady state (0.3 mM/s, from Ref. 21; and 0.4 mM/s, from Ref. 6). The estimate for the Hill coefficient, *n*H, was on average 2.57 ± 0.17. The fact that the final estimate is not very different from the original prior suggests that the Hill function adequately explains the dynamics of the system during the experiment. The estimate for *K*
_{ADP} is 58 μM, similar to (30% larger than) the value reported previously (21).

The cellular ATPase rate, from rest to activity, increased ∼10-fold for these moderate stimulation rates (1 Hz), which in this muscle reduced the PCr concentration to about one-half of its resting value. From the basal ATPase rate (0.0117 mM/s) and the value for*V*
_{max} (0.50 mM/s), the maximal change of OxPhos (i.e., basal-to-maximal flux) can be estimated at ∼43-fold. This value is somewhat less than estimated from changes in blood flow in human quadriceps musculature (2), but the forearm musculature in the untrained, largely sedentary subjects used here contains fewer oxidative cells than in the thigh, so a lower value is expected. These values confirm quantitatively the qualitative interpretation made in Ref. 6, that the muscles operated within their capacity to attain energy balance. That is, the increase in ATPase required for the 1-Hz twitch stimulation used in the original experiments, on average, did not metabolically stress the energetic system to its limit. Finally, the average coefficient of variations for individual estimates of the parameters were always much smaller than the percent dispersion of the average of the group. This result confirms the earlier result that there are significant inter-individual differences in the metabolic characteristics relevant to energy balance in normal subjects (7).

The conclusion from this phase of the analysis is that the basic*Model 1* in *Eqs. 3-6
* accounts reasonably well for energy balance; that is, the model fits the trend of the data from this somewhat complex experimental design, suggesting that the three components given in *Eq. 3
* and their particular kinetic expressions are necessary to account for the integrated operation of intracellular energy metabolism. Moreover, the biochemical parameters are in the range of expected values given in the literature.

However, despite the overall quality of the model fits, Fig. 3 reveals systematic deviations and residuals that demand further exploration. The next sections test two hypotheses about the sufficiency of the model and whether inclusion of specific additional features into the basic energy balance model can reduce these deviations.

### Model 2: Role of Cell Type Heterogeneity

The model residuals (i.e., the difference between model prediction and PCr and P_{i} time courses) were not normally distributed during recovery. One explanation is a variation in the mitochondria enzyme activity and percentage volume of mitochondria in different fibers; the initial phase of the recovery would be dominated by those cells with more mitochondria. Alternatively, the ATPase fluxes during stimulation may be different in different fibers, because of differences in actomyosin and sarcoplasmic reticular ATPase isoforms, to list the two largest ATPases in the cell. We now test the hypothesis of whether heterogeneity of cellular properties might reduce the differences between model prediction and data. Heterogeneity of mitochondrial capacity among the cells in the forearm musculature has been reported in some individuals (34) and is consistent with analyses of cell types by histochemical staining of the forearm wrist and finger flexors (22).

The forearm flexor muscle studied in the original experiments consisted of a range of muscle cell types among the different individuals (7), but the importance of heterogeneity within an individual was not assessed. The significant question addressed now is whether heterogeneous properties at the cellular level are sufficient to make a difference in the macroscopic observations of the composite muscle. That differing fiber types may compose muscles is well known (35). Human muscle cells have a range of mechanical speeds (28). However, histochemically identified type 1 and type 2 fibers have no difference in concentration of relevant metabolites (18), in contrast with the large differences in animal muscle (27).

We analyzed the possibility that accounting for differences in type 1 and type 2 fibers' energetics would improve model performance by defining two parallel pathways in the model, each with different values of *k*
_{stim} for contractile ATPase and*V*
_{max} for OxPhos flux. The fraction of slow type 1 fibers was adjustable within Bayesian constraints (0.47 ± 0.10) using available prior knowledge (22). This implementation of the results in Ref. 22 allows the fitting procedure to adjust individually the type 1 fiber fraction, based on the individual subject's data. Stienen et al. (42) showed the relative difference in ATPase between slow and fast fibers was ∼2.5-fold, and this factor was used to constrain the model:*k*
_{stim}
^{f} = 2.5 ·*k*
_{stim}
^{s}, where the superscripts “f” and “s” indicate slow and fast fibers, respectively. There is an approximate twofold difference in the content of mitochondria between the more and less oxidative cells, so this factor was also used in the model: *V*
_{max}
^{s} = 2 · *V*
_{max}
^{f}. Table5 gives the salient parameter values thus found. The fraction of fast fibers found a posteriori in the fitting averaged 48% (with individual values ranging from 24 to 61%), which is not different statistically from that reported in Ref. 22.

The results of this analysis shows that the value for*K*
_{ADP} becomes substantially lower (from 58 ± 21 to 36 ± 10 μM) and the value for *n*H substantially higher (from 2.57 ± 0.17 to 3.43 ± 0.78) than with the basic model. Our model for heterogeneity assumed that mitochondria in both types were qualitatively identical; only the number (and thus the activity of OxPhos flux enzymes) differed. Estimated *k*
_{rest} remains somewhat higher than expected, and the *V*
_{max} values for oxidative phosphorylation and *k*
_{stim} are within expectation. It is of note that, averaging the values of slow and fast fibers, one obtains 0.495 mM/s and 0.0142 1/s for*V*
_{max} and *k*
_{stim}respectively, which are not different from the values found with the homogeneous *Model 1* and are consistent with an average 50% fiber distribution. Comparison of Figs. 3 and4 demonstrates some improvement in the quality of fit during the resting phase of the experiment but similar deviations during recovery. AIC and BIC values show (see Table 7) that this model accounts better for the data than *Model 1*. These model results suggest that fiber type distribution is probably among the sources of the considerable inter-individual differences in parameter estimates reported earlier for these same subjects (7) and confirmed here.

### Model 3: Analysis of a Different Control of Metabolic Rate During Resting and Contractile States

In this section, we propose and test the concept that the control of oxidative phosphorylation is quantitatively different at rest from the control during activity and recovery. This hypothesis suggested itself because model estimates of basal ATPase flux are always higher than expected. The basal ATPase flux estimated in Table 4 (0.012 mM/s) is almost twofold higher than measured in the analyses cited in Table 3(0.007 mM/s). The model predictions and data also consistently differed during the initial resting phase (Figs. 3 and 4). Thus these results raise the possibility that *Models 1* and *2* may not apply equally well to muscle at rest and during activity.

The basic premise of our model development is that ADP is the molecular signal, that is, shows feedback regulation, and, like other cellular signals, it has characteristics of ultrasensitivity (23, 24). Much evidence points to the possibility of a feed-forward type of regulation operating in parallel, specifically involving Ca^{2+} activation of mitochondrial enzymes (see Refs. 16 and 17 for a review). We tested this class of mechanism with our model by including a new parameter, φ, as a multiplier of *V*
_{max} in the OxPhos flux expression (*Eq. 4
*). The magnitude of the multiplier depends on the state of the muscle: its value is unitary for unstimulated resting state and higher during activity and recovery phases (thus accounting for the possibility that *V*
_{max} increased in response to contractile activity by a feed-forward mechanism). This resulted in modification of *Eq. 4
*
Equation 15We chose φ = 3 after preliminary simulations where φ varied between 1 and 10. In the modified model, this feed-forward activation (increase in the value) of *V*
_{max} is added during the recovery phase (*phase E*) only, and not at rest (during *phases A*, *B*, and*F*). Parameter results are shown in Table6. The quality of fit improved substantially (Fig. 5), primarily during the resting phases of the experiment (often leading to lower estimates of *k*
_{rest} close to the independently measured value) and during the recovery phase. The model resulted in slightly lower estimates for *K*
_{ADP} without major differences in *k*
_{stim} (contractile ATPase flux). Overall, the parameter values were similar to those found in the other two models. We also tested feed-forward increases in*K*
_{ADP} with similar results; however,*K*
_{ADP} was more sensitive to the magnitude of φ than was *V*
_{max}, because the term for*K*
_{ADP} has the exponent *n*H. AIC and BIC analysis (Tables 7 and 8) suggests that*Model 3* is superior both to *Model 1* and*Model 2* in terms of data descriptive power and parsimony.

This model analysis suggests that properly designed experiments can detect differences in *V*
_{max} of OxPhos flux between resting and active states: a combined feedback and feed-forward control of OxPhos flux may be necessary to account for the energetics of skeletal muscle at rest and during activity, with the same parameters in one model. Thus it may not be valid to assume that the regulation of OxPhos flux occurs by the same control mechanisms over the entire range of cellular respiration.

### Role of Glycolysis

Recent work demonstrates that glycogenolysis and glycolysis contribute to energy balance in the human limb muscle (13,14), but there are no terms for glycolysis in our model. Glycolytic ATP synthesis is implicitly included in the other terms for ATP synthesis (creatine kinase reaction transiently and net synthesis by oxidative phosphorylation). This is one reason we know the present models are incomplete. For example, the model-estimated ATPase tended toward zero during *phase D* (Fig. 1), whereas it was expected to be similar to the basal state, if not higher, due to ion-pumping energy requirements immediately after the stimulation. Significant and unaccounted ATP synthesis during *phase D* (Fig. 1), such as glycolysis, could explain this result.

It is easy to show that these models cannot account for the pH changes reported in the original reports (6, 13,14). In the experiments modeled, the pH response is biphasic, with an initial alkalization during the initial phase of muscle contraction followed by an acidification that lasts well into recovery. Simulations with the model reactions show a monotonic alkalization with PCr decrease and its reversal when PCr is restored. Thus there is a marked difference between the measured pH and that predicted from the PCr(*t*) calculated by *Model 1*(data not shown). These differences are due to lactate production in glycolytic ATP synthesis (13, 14). These comparisons demonstrate that our models are sensitive to the absence of a term for glycolysis and likely could better account for the data if*Eq. 5
* included appropriate terms. We have not yet accomplished this task, because a kinetic expression for the control of glycolysis and glycogenolysis is not yet available. High-quality NMR data are needed analyze the time course of pH changes, because the pH changes are quite small in these relatively low-intensity exercises.

## DISCUSSION

A prerequisite to understanding the regulation of cellular energy balance for muscle or any other cell is developing a mechanistic and quantitative description of the processes involved. The simple models reported here use a few kinetic mechanisms established in the literature, some from in vitro and others from in vivo experiments, united in a conceptually simple set of ordinary differential equations. Despite the highly nonlinear properties of these equations and the low number of mechanisms modeled, the resulting final picture accurately accounts for the major features of a complex and energetically significant muscle stress, both at rest and during activity, involving both ischemia and normal blood flow. We know that significant components of bioenergetics, such as glycolytic ATP synthesis and Ca^{2+} influences on mitochondrial enzyme activity, are missing. The results of the model analyses reported here would likely be changed were these mechanisms included. In this sense, the models we describe here are wrong and incomplete. However, the concepts and mechanisms involved are well established, and in fact, we build on early advances. Our analysis of energy balance is conceptually similar to the energy-balance model proposed by Meyer and Foley (31, 33). Our model contains kinetic expressions for the forward and reverse kinetics of this reaction, and the net flux is explicitly the difference between the two unidirectional fluxes that are determined by the time-dependent metabolite concentration. Furthermore, the form of Meyer's analysis yields a mono-exponential time course for PCr decline and recovery. Our present model does not so constrain the time course; nonetheless, the properties of the creatine kinase reaction lead to a time course very close to a single exponential (25). Similarly, we consider ADP concentration as the control for oxidative phosphorylation, as did Ref. 9. However, we used a different kinetic expression, given in Ref.21. Thus we emphasize that our current modeling approach does allow tests of current knowledge of the system, because they only contain the parsimonious representation of functions thought to be correct and dominant in energy balance. In this sense, our work establishes new insights and criteria for analyzing experimental data. We make explicit and mechanistic formulations of the processes involved. These can readily be changed as additional information is made available. For example, the ATPase rate in the current model is a continuous steady rate during contractile activity; intermittent maximal electrical twitches are often used and were used in the experiments analyzed here, making the ATPase episodic. We have not yet fully analyzed the effect this characteristic would have; however, the effect is expected to be small because of the widely different characteristic times of the processes involved: less than 1 s for the mechanically driven ATPases and many seconds for the mitochondria ATP synthesis. Thus the metabolic response time may be too slow to reflect the rapid transients in the ATPase. In addition, it is possible to measure the force production or work done by the muscle, and changes in this output could be factored into the expression for contractile ATPase, for example, testing the hypothesis that the chemical cost for a given biomechanical output is independent of the chemical changes occurring in the cell.

An important aspect of our work is the explicit use of statistical tools that provide a method to assess whether the inclusion of different features and components of the model renders a better description of the observed data. The AIC and BIC criteria shown in Table 7 are to our knowledge the best for this purpose. These criteria demonstrate that our third model best describes the data from these experiments. The improvement in explanatory power of a model derived from adding a certain mechanistic component can be quantified statistically based upon the criteria that we employed. Differences in AIC and BIC can then be appropriately tested for significance (Table 8), thus giving a basis for hypothesis testing. The indication would be that either the added complexity is unwarranted or that it results in a significant improvement in the model behavior. Further work will indicate whether this conclusion is valid for other data sets and whether inclusion of other components discussed above will actually improve the fit. These criteria also show limits to making experimental conclusions, because, if the model is conceptually correct, as we believe is the case, then addition of a putative “essential” mechanism should improve the fit. If that addition has independently verifiable validity, but the model fit to the data is not improved, then it follows that a better experimental design and/or precision is needed to demonstrate the effects of that component in the intact system.

There are several physiological implications in our results; however, as discussed in the next paragraph, we cannot yet take these results as established. *Model 1* considers that all cells are identical in their biochemical properties and that the same function for ADP control of oxidative phosphorylation applies to resting and contracting muscle. The fit is rather good, considering the simplicity of the model and the potential complexity of the intracellular environment (localization of relevant enzymes and potential compartmentalization of metabolites and reactions). All things considered, metabolites freely diffusing to enzymes is a reasonable starting point for quantitative analysis of intracellular energy metabolism. This same assumption is made in the other two models, and our results are consistent with recent evidence from radial and longitudinal diffusivity of relevant metabolites in muscle (19, 20). *Model 2* explicitly considers the two cell types normally found in mammalian muscle: fast and slow. The fit to the data showed only a small, but significant, improvement (Tables7 and 8). Thus we can conclude that analysis of muscle as a single cell type is useful as an approximation but in detail is insufficient; however, we note that the degree of insufficiency as judged from Figs. 3 and 4 is small and may be difficult to detect. Note that our analysis is the first to unite the resting state with the contractile state; in classic energetic analyses, the contractile state is always superimposed on a resting baseline which is subtracted and not quantitatively analyzed. A bigger improvement was seen with *Model 3*, in which we hypothesized that the control of oxidative phosphorylation during contractile activity was qualitatively different than at rest. Thus analysis of the data with *Model 3* indicated the flux of oxidative phosphorylation might have dual control. Control can largely be attributed to the ultrasensitive signal function of ADP (21) or to related energetic functions as the chemical potential of ATP (33) as the results with *Models 1* and *2* show, but this control appears to be insufficient to account for both the dynamic changes in PCr and P_{i} and their resting values. A feed-forward factor substantially improves the fit (Tables 7 and 8). This modeling result does not, of course, establish a dual control mechanism, but this novel result has wide implications that warrant further experimental analyses.

Model validation is a complex issue, and it entails the careful comparison of model-estimated quantities with other, independently determined, measures of the same quantities. We have not attempted validation here, because it is impossible with these data alone. However, the results of our analysis suggest directions for future experimental design that would help validating such a modeling construct and interpreting the parameter estimates as physiologically relevant parameters. It is obviously a challenge to obtain sufficient human muscle material to make the appropriate characterizations of the mitochondria parameters (*K*
_{ADP},*V*
_{max}, *n*H) in vitro. Additional data in the literature when analyzed, including animal experiments, would test whether a consistent set of biochemical parameters that had been measured independently would be obtained. Such results would provide a measure of validation. Human subjects can be selected according to physical activity, and animals can be similarly selected after known adaptive procedures, and both can be biopsied for independent analyses. Models such as this would then be potentially useful to provide estimates of important physiological and biochemical properties of individual subjects, since the basic NMR acquisition is noninvasive. This approach would make individualized bioenergetic parameter estimates feasible. Knowledge of important physiological and biochemical properties may find clinical applications in documenting the progression of certain chronic disease states of muscle, judging the utility of therapeutic interventions, and assessing the success of adaptive strategies in sports medicine. These opportunities provide strong motivation for the continued experimental and model evaluations of cellular energy metabolism.

Another limitation in the data obtained from intact muscle by NMR spectroscopy is that not all parameters are separated experimentally equally well, i.e., there are distinct phases in the experiment where individual kinetic processes dominate with respect to others: for example, during ischemic stimulation, contractile ATPase dominates and oxidative phosphorylation is absent. During recovery, oxidative phosphorylation dominates when ATPase reverts to basal. As expected, model fitting results show no correlation between these two parameters. On the other hand, the terms for *K*
_{ADP},*V*
_{max}, and *n*H appear in the same kinetic expression for oxidative phosphorylation, and, as would be expected, the model analyses show considerable correlation between the magnitude of these terms: *K*
_{ADP} and*V*
_{max} are always positively correlated (correlation coefficients vary between 0.7 and 0.9 in different individuals and models), and *V*
_{max} and*n*H are always negatively correlated (correlation coefficients vary between −0.4 and −0.6 in the different individuals and models). It must be noted that high correlations, however, do not always result in a more difficult model identification. From our results, it seems that all model parameters are well resolved (or, equivalently, that the model prediction is reasonably sensitive to their value). Thus more complicated experiments that would separately manipulate those parameters, such as through known genetic states, would be valuable.

The types of data available by ^{31}P NMR methods also drive our “minimal modeling” philosophy. The present model accommodates only those mechanisms that are likely to affect the data. Thus the resulting estimates of kinetic parameters are “aggregate” or “compound” quantities, which implicitly include the effect of other processes going on at the same time. For example, we set ATPase to zero during the poststimulation ischemic interval (*phase D*, Fig.1), yet, it is likely that glycolysis occurs during this phase (see*Role of Glycolysis*, above in results). Thus the apparent absence of an ATPase term during this phase means that two competing processes, not specified in the model, occur simultaneously. Such components can be added to the model, provided additional experimental data reflecting those mechanisms are available (analysis of muscle ionic content for the former and of lactate for the latter) and appropriate mechanistic equations are derived.

The weakness of our simplified modeling of macroscopic observations is that the details or even the existence of additional mechanisms cannot always be interpreted directly, but other types of investigation have done that. The main strength of our approach is its integration of known mechanisms and the ease of extending it by inclusion of additional components. For example, we saw that we could analyze an additional piece of data, like pH, in the context of the model and deduced the need for an additional glycolytic ATP synthesis. The great success of past decades of mechanistic and often in vitro analyses of components, well-defined in molecular terms, gives us a large compendium of specific mechanisms to consider for the analysis of energy balance. The fact that our choice of mechanisms accounts for so much of the observed quantities means that we have identified the quantitatively most significant components. Our “simple model” is therefore quite complex in the sense usually applied to living systems; several simple equations or “rules” explain highly complex behavior.

The present model analysis also demonstrates several ways by which data collection might improve. Higher time resolution will always be helpful in NMR spectroscopy. The timing of the physiological events needs to be as accurate as possible. In the experiments analyzed here (6), the timing was noted only to the interval within one spectrum, ±9 s; the timing error would be significantly greater with any miscounting of spectral acquisition. Our analysis with this model shows that the parameter estimates are highly sensitive to event timing, e.g., when increase in muscle ATPase occurs during activity or when restoration of blood flow occurs. The interaction between the accuracy of peak integration from the NMR spectra and the number of metabolite concentration points per time interval has not been investigated; in general, the shorter the time interval, the greater the variance in peak area estimates. It may prove desirable to acquire high time resolution only during the crucial times of the experiment, when physiological events change rapidly, and acquire spectra with a higher signal-to-noise ratio at other times.

## Acknowledgments

We gratefully acknowledge the editorial help of Eileen Thorsos.

## Footnotes

This work was partially supported by National Institutes of Health Grants NCRR-12609, AR-41928, and AR-36281.

Address for reprint requests and other correspondence: P. Vicini, Dept. of Bioengineering, Box 352255, Univ. of Washington, Seattle, WA 98195-2255 (E-mail: vicini{at}u.washington.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

- Copyright © 2000 the American Physiological Society