Departments of 1 Bioengineering, 2 Radiology, and 3 Physiology and Biophysics, University of Washington, Seattle, Washington 98195
![]() |
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 (Pi) data obtained in normal subjects with an activity-recovery protocol using 31P 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 Pi 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
![]() |
INTRODUCTION |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
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. 31P 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
![]() |
(1) |
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 Pi. 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. Table 1 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 31P 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 (see DISCUSSION).
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 Pi. 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 Pi (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 (ADP0) and ATP
(ATP0) concentrations and the ATP concentration at
time t
![]() |
(2) |
The following equation describes the change in ATP concentration over
time
![]() |
(3) |
A Hill sigmoid function describes the time course of the flux of ATP
synthesis by OxPhos flux (21)
![]() |
(4) |
The equation describing PCr change over time is
![]() |
(5) |
The following equation describes the time course of Pi
during PCr and ATP transient
![]() |
(6) |
Parameter Estimation
Predicted PCr concentration, PCr(t), the solution to Eq. 5, and Pi concentration, Pi(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)
![]() |
(7) |
![]() |
(8) |
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
![]() |
(9) |
After defining the MAP objective function, we performed parameter
estimation on the entire set of parameters Np.
Note that it is still possible to calculate the precision of the
parameter estimates using an extension of Eq. 8
![]() |
(10) |
Experimental evidence suggested that the measurement errors affecting the Pi and PCr data were assumed to be Gaussian, zero mean, and with a constant standard deviation (Pi) 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 Pi peak integrals increases as Pi increases. Thus the measurement error will change during the different phases of the experiment as PCr and Pi 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 Table 2 are well known from a priori information. ATP0 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), PCr0, and the value for ADP0 can both be uniquely calculated as a function of other model parameters from the steady-state equations
![]() |
(11) |
![]() |
(12) |
![]() |
(13) |
![]() |
(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): KADP, nH, and krest. Prior mean and standard deviation values were taken from Ref. 21 for KADP = 0.043 ± 0.003 mM and for nH = 2.11 ± 0.34. Mean and standard deviation for krest (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. Vmax and kstim 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. Table
3 displays the numerical results of this
analysis. These values for krest yield an
average basal ATPase of 0.007 ± 0.001 mM/s (0.00087 s1 · 8.2 mM), consistent with the result given in
Ref. 6, cited above. This analysis of basal ATPase flux can constrain
krest to follow a certain mean and standard
deviation. Therefore, we used a Bayesian prior on
krest 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, PCr0, and the initial ADP concentration, ADP0, were all in the expected physiological range. On average, these were 0.016 ± 0.008 mM for ADP0 and 32.1 ± 3.7 mM for PCr0. The estimated inorganic phosphate initial concentration, Pi0, 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 Vmax from these transient
experiments is 0.50 ± 0.34 mM ATP/s, somewhat larger than the
estimates for Vmax 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, nH, 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 KADP 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 Vmax (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 Pi 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 kstim for contractile ATPase and
Vmax 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: kstimf = 2.5 · kstims, 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: Vmaxs = 2 · Vmaxf. Table
5 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
KADP becomes substantially lower (from 58 ± 21 to 36 ± 10 µM) and the value for nH
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 krest remains somewhat higher than
expected, and the Vmax values for oxidative
phosphorylation and kstim 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 Vmax and kstim
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 and
4 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 Ca2+ 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 Vmax 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 Vmax increased in
response to contractile activity by a feed-forward mechanism). This
resulted in modification of Eq. 4
![]() |
(15) |
|
|
|
This model analysis suggests that properly designed experiments can detect differences in Vmax 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 Ca2+ 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 (Tables
7 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
Pi 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 (KADP, Vmax, nH) 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 KADP,
Vmax, and nH 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: KADP and
Vmax are always positively correlated
(correlation coefficients vary between 0.7 and 0.9 in different
individuals and models), and Vmax and
nH 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 31P 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.
![]() |
ACKNOWLEDGEMENTS |
---|
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.
Received 19 April 1999; accepted in final form 4 January 2000.
![]() |
REFERENCES |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
1.
Akaike, H.
Information theory and an extension of the maximum likelihood principle.
In: 2nd International Symposium on Information Theory, edited by Petrov BN,
and Caski F.. Tsahkadsor, Armenia: USSR Acad Sci, 1971, p. 267-281.
2.
Anderson, P,
and
Saltin B.
Maximal perfusion of skeletal muscle in man.
J Physiol (Lond)
366:
233-249,
1985[Abstract].
3.
Awad, AM.
Properties of the Akaike Information Criterion.
Microelectronics Reliability
36:
457-464,
1996.
4.
Barrett, PH,
Bell BM,
Cobelli C,
Golde H,
Schumitzky A,
Vicini P,
and
Foster DM.
SAAM II: simulation, analysis, and modeling software for tracer and pharmacokinetic studies.
Metabolism
47:
484-492,
1998[ISI][Medline].
5.
Bell, BM,
Burke JV,
and
Schumitzky A.
A relative weighting method for estimating parameters and variances in multiple data sets.
Comp Stats Data Anal
22:
119-135,
1996.
6.
Blei, ML,
Conley KE,
and
Kushmerick MJ.
Separate measures of ATP utilization and recovery in human skeletal muscle.
J Physiol (Lond)
465:
203-222,
1993[Abstract].
7.
Blei, ML,
Conley KE,
Odderson IB,
Esselman PC,
and
Kushmerick MJ.
Individual variation in contractile cost and recovery in a human skeletal muscle.
Proc Natl Acad Sci USA
90:
7396-7400,
1993[Abstract].
8.
Carson, ER,
Cobelli C,
and
Finkelstein L.
The Mathematical Modeling of Metabolic and Endocrine Systems. New York: Wiley, 1982.
9.
Chance, B,
Leigh JS, Jr,
Kent J,
and
McCully K.
Metabolic control principles and 31P NMR.
Fed Proc
45:
2915-2920,
1986[ISI][Medline].
10.
Chance, B,
Leigh JS,
Kent J,
McCully K,
Nioka S,
Clark BJ,
and
Maris JM.
Multiple controls of oxidative metabolism in living tissues as studied by phosphorus magnetic resonance.
Proc Natl Acad Sci USA
83:
9458-9462,
1986[Abstract].
11.
Cobelli, A,
Caumo A,
and
Omenetto M.
Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian two-compartment model.
Am J Physiol Endocrinol Metab
277:
E481-E488,
1999
12.
Cobelli, C,
and
DiStefano JJ.
Parameter and structural identifiability concepts and ambiguities: a critical review and analysis.
Am J Physiol Regulatory Integrative Comp Physiol
239:
R7-R24,
1980[ISI][Medline].
13.
Conley, KE,
Blei ML,
Richards TL,
Kushmerick MJ,
and
Jubrias SA.
Activation of glycolysis in human muscle in vivo.
Am J Physiol Cell Physiol
273:
C306-C315,
1997
14.
Conley, KE,
Kushmerick MJ,
and
Jubrias SA.
Glycolysis is independent of oxygenation state in stimulated human skeletal muscle in vivo.
J Physiol (Lond)
511:
935-945,
1998
15.
Funk, CI,
Clark AJ,
and
Connett RJ.
A simple model of aerobic metabolism: applications to work transitions in muscle.
Am J Physiol Cell Physiol
258:
C995-C1005,
1990
16.
Gunter, TE,
Gunter KK,
Sheu SS,
and
Gavin CE.
Mitochondrial calcium transport: physiological and pathological relevance.
Am J Physiol Cell Physiol
267:
C313-C339,
1994
17.
Hansford, RG.
Role of calcium in respiratory control.
Med Sci Sports Exerc
26:
44-51,
1994[ISI][Medline].
18.
Harris, RC,
Hultman E,
and
Nordesjo LO.
Glycogen, glycolytic intermediates and high-energy phosphates determined in biopsy samples of musculus quadriceps femoris of man at rest. Methods and variance of values.
Scand J Clin Lab Invest
33:
109-120,
1974[ISI][Medline].
19.
Hubley, MJ,
Locke BR,
and
Moerland TS.
Reaction-diffusion analysis of the effects of temperature on high-energy phosphate dynamics in goldfish skeletal muscle.
J Exp Biol
200:
975-988,
1997
20.
Hubley, MJ,
Rosanske RC,
and
Moerland TS.
Diffusion coefficients of ATP and creatine phosphate in isolated muscle: pulsed gradient 31P NMR of small biological samples.
NMR Biomed
8:
72-78,
1995[ISI][Medline].
21.
Jeneson, JAL,
Wiseman RW,
Westerhoff HV,
and
Kushmerick MJ.
The signal transduction function for oxidative phosphorylation is at least second order in ADP.
J Biol Chem
271:
27995-27998,
1996
22.
Johnson, MA,
Polgar J,
Weightman D,
and
Appleton D.
Data on the distribution of fibre types in thirty-six human muscles: an autopsy study.
J Neurol Sci
18:
111-129,
1973[ISI][Medline].
23.
Koshland, DE, Jr.
Switches, thresholds and ultrasensitivity.
Trends Biochem Sci
12:
225-229,
1987[ISI].
24.
Koshland, DE, Jr,
Goldbeter A,
and
Stock JB.
Amplification and adaptation in regulatory and sensory systems.
Science
217:
220-225,
1982[ISI][Medline].
25.
Kushmerick, MJ.
Energy balance in muscle activity: simulations of ATPase coupled to oxidative phosphorylation and to creatine kinase.
Comp Biochem Physiol B Biochem Mol Biol
120:
109-123,
1998[ISI][Medline].
26.
Kushmerick, MJ.
Multiple equilibria of cations with metabolites in muscle bioenergetics.
Am J Physiol Cell Physiol
272:
C1739-C1747,
1997
27.
Kushmerick, MJ,
Moerland TS,
and
Wiseman RW.
Mammalian skeletal muscle fibers distinguished by contents of phosphocreatine, ATP, and Pi.
Proc Natl Acad Sci USA
89:
7521-7525,
1992[Abstract].
28.
Larsson, L,
and
Moss RL.
Maximum velocity of shortening in relation to myosin isoform composition in single fibres from human skeletal muscles.
J Physiol (Lond)
472:
595-614,
1993[Abstract].
29.
Lawson, JWR,
and
Veech RL.
Effects of pH and free Mg2+ on the Keq of the creatine kinase reaction and other phosphate hydrolysis and phosphate transfer reactions.
J Biol Chem
254:
6528-6537,
1979[Abstract].
30.
McFarland, EW,
Kushmerick MJ,
and
Moerland TS.
Activity of creatine kinase in a contracting mammalian muscle of uniform fiber type.
Biophys J
67:
1912-1924,
1994[Abstract].
31.
Meyer, RA.
A linear model of muscle respiration explains monoexponential phosphocreatine changes.
Am J Physiol Cell Physiol
254:
C548-C553,
1988
33.
Meyer, RA,
and
Foley JM.
Cellular processes integrating the metabolic response to exercise.
In: Handbook of Physiology. Exercise: Regulation and Integration of Multiple Systems. Bethesda, MD: Am. Physiol. Soc, 1996, sect. 12, part III, chapt. 18, p. 841-869.
34.
Mizuno, M,
Secher NH,
and
Quistorff B.
31P-NMR spectroscopy, rsEMG, and histochemical fiber types of human wrist flexor muscles.
J Appl Physiol
76:
531-538,
1994
35.
Pette, D,
and
Staron RS.
The molecular diversity of mammalian muscle fibers.
News Physiol Sci
8:
153-157,
1993
36.
Radda, GK.
The use of NMR spectroscopy for the understanding of disease.
Science
233:
640-645,
1986[ISI][Medline].
37.
Sakamoto, Y,
Ishiguro I,
and
Kitagawa G.
Akaike Information Criterion Statistics. Tokyo, Japan: KTK Scientific, 1986.
38.
Schimerlik, ML,
and
Cleland WW.
Inhibition of creatine kinase by chromium nucleotides.
J Biol Chem
248:
8418-8423,
1973
39.
Schwarz, G.
Estimating the dimension of a model.
Annu Stat
6:
461-464,
1978.
40.
Sjogaard, G,
and
Saltin B.
Extra- and intracellular water spaces in muscles of man at rest and with dynamic exercise.
Am J Physiol Regulatory Integrative Comp Physiol
243:
R271-R280,
1982
41.
Smeitink, J,
Wevers R,
Hulshof J,
Ruitenbeek W,
Lith TV,
Sengers R,
Trijbels F,
Korenke C,
and
Walliman T.
A method for quantitative measurement of mitochondrial creatine kinase in human skeletal muscle.
Ann Clin Biochem
20:
196-201,
1992.
42.
Stienen, GJM,
Kiers JL,
Bottinelli R,
and
Reggiani C.
Myofibrillar ATPase activity in skinned human skeletal muscle fibres: fibre type and temperature dependence.
J Physiol (Lond)
493:
299-307,
1996[Abstract].
43.
Zurlo, F,
Nemeth PM,
Choksi RM,
Sesodia S,
and
Ravussin E.
Whole-body energy metabolism and skeletal muscle biochemical characteristics.
Metabolism
43:
481-486,
1994[ISI][Medline].