Cellular energetics analysis by a mathematical model of energy balance: estimation of parameters in human skeletal muscle

Paolo Vicini1 and Martin J. Kushmerick1,2,3

Departments of 1 Bioengineering, 2 Radiology, and 3 Physiology and Biophysics, University of Washington, Seattle, Washington 98195


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
REFERENCES

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
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
REFERENCES

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
PCr<SUP><IT>−</IT></SUP><IT>+</IT>MgADP<SUP><IT>−</IT></SUP><IT>+</IT>H<SUP><IT>+</IT></SUP><IT> ↔ </IT>MgATP<SUP><IT>−</IT></SUP><IT>+</IT>Cr (1)
The 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 Pi, 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 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
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
REFERENCES

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.

                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Summary of the changes of model parameters during the different phases of the experiment design schematically diagrammed in Fig. 1



View larger version (12K):
[in this window]
[in a new window]
 
Fig. 1.   Schematic of experimental (QUEST) protocol. See text for details. Phases A-F are indicated.

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 Delta 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
&Dgr;ADP(<IT>t</IT>)<IT>=</IT>ADP<SUB><IT>0</IT></SUB><IT>+</IT>ATP<SUB><IT>0</IT></SUB><IT>−</IT>ATP(<IT>t</IT>) (2)
Three 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
<FR><NU>dATP(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=</IT>−<IT>k</IT>ATP(<IT>t</IT>)<IT>+</IT>OxPhos(<IT>t</IT>)<IT>−</IT><FR><NU>dPCr(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR> (3)
where PCr0 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)/dt 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, krest, to a different stimulated value, kstim, 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)
OxPhos(<IT>t</IT>)<IT>=</IT><FR><NU><IT>V</IT><SUB>max</SUB><FENCE><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ADP</SUB></DE></FR></FENCE><SUP><IT>n</IT>H</SUP></NU><DE><IT>1+</IT><FENCE><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ADP</SUB></DE></FR></FENCE><SUP><IT>n</IT>H</SUP></DE></FR> (4)
where Vmax is the maximal velocity of OxPhos flux at infinite ADP concentration, KADP is the substrate concentration at which the flux of ATP synthesis by OxPhos flux equals half of Vmax and nH is the Hill coefficient.

The equation describing PCr change over time is
<FR><NU>dPCr(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><FR><NU>−<FR><NU>[<IT>&Dgr;</IT>ADP(<IT>t</IT>)]PCr(<IT>t</IT>)<IT>V</IT><SUB>for CK</SUB></NU><DE><IT>K</IT><SUB>b</SUB><IT>K</IT><SUB>ia</SUB></DE></FR><IT>+</IT><FR><NU>ATP(<IT>t</IT>)[TCr<IT>−</IT>PCr(<IT>t</IT>)]<IT>V</IT><SUB>rev CK</SUB></NU><DE><IT>K</IT><SUB>iq</SUB><IT>K</IT><SUB>p</SUB></DE></FR></NU><DE><IT>1+</IT><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ia</SUB></DE></FR><IT>+</IT><FR><NU>ATP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>iq</SUB></DE></FR><IT>+</IT><FR><NU>PCr(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ib</SUB></DE></FR><IT>+</IT><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)PCR(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>b</SUB><IT>K</IT><SUB>ia</SUB></DE></FR><IT>+</IT><FR><NU>ATP(<IT>t</IT>)[TCr<IT>−</IT>PCr(<IT>t</IT>)]</NU><DE><IT>K</IT><SUB>iq</SUB><IT>K</IT><SUB>p</SUB></DE></FR></DE></FR> (5)
The constants Kia, Kiq, Kib, Kb, and Kp 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. Vfor CK is the maximal velocity in the forward CK flux, and Vrev 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 Pi during PCr and ATP transient
P<SUB>i</SUB>(<IT>t</IT>)<IT>=</IT>P<SUB>i<IT>0</IT></SUB><IT>+</IT>[ATP<SUB><IT>0</IT></SUB><IT>−</IT>ATP(<IT>t</IT>)]<IT>+</IT>[PCr<SUB><IT>0</IT></SUB><IT>−</IT>PCr(<IT>t</IT>)] (6)
where Pi0 is the steady-state value at rest for Pi. 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 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)
WRSS(<A><AC>p</AC><AC>ˆ</AC></A>)<IT>=</IT><LIM><OP>∑</OP><LL><IT>i=1</IT></LL><UL><IT>N</IT></UL></LIM> <FR><NU>[<IT>y<SUB>i</SUB>−</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t<SUB>i</SUB></IT>)]<SUP><IT>2</IT></SUP></NU><DE>V<SUB><IT>i</IT></SUB></DE></FR> (7)
where N is the number of data points, yi is the data point at time ti, p is the estimated parameter vector (size Np), m(p,ti) is the model function evaluated at the parameter vector p and ti, and Vi is the measurement error variance of the data point yi. 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
Cov(<A><AC>p</AC><AC>ˆ</AC></A>)<IT>≥</IT><FENCE><FENCE><FR><NU><IT>∂</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t</IT>)</NU><DE><IT>∂</IT><A><AC>p</AC><AC>ˆ</AC></A></DE></FR></FENCE><SUP>T</SUP>diag<FENCE><FR><NU><IT>1</IT></NU><DE>V<SUB><IT>i</IT></SUB></DE></FR></FENCE> <FR><NU><IT>∂</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t</IT>)</NU><DE><IT>∂</IT><A><AC>p</AC><AC>ˆ</AC></A></DE></FR></FENCE><SUP><IT>−1</IT></SUP> (8)
where 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
MAP(<A><AC>p</AC><AC>ˆ</AC></A>)<IT>=</IT><LIM><OP>∑</OP><LL><IT>i=1</IT></LL><UL><IT>N</IT></UL></LIM> <FR><NU>[<IT>y<SUB>i</SUB>−</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t<SUB>i</SUB></IT>)]<SUP><IT>2</IT></SUP></NU><DE>V<SUB><IT>i</IT></SUB></DE></FR><IT>+</IT><LIM><OP>∑</OP><LL><IT>i=1</IT></LL><UL><IT>N</IT><SUB>b</SUB></UL></LIM> <FR><NU>[<IT>&mgr;<SUB>i</SUB>−</IT><A><AC>p</AC><AC>ˆ</AC></A><SUB><IT>i</IT></SUB>]<SUP><IT>2</IT></SUP></NU><DE><IT>&Sgr;<SUB>i</SUB></IT></DE></FR> (9)
where Nb is the number of (uncorrelated) Bayesian parameters (smaller than or equal to the total number of parameters in the model Np), µi is the mean, Sigma i is the variance associated with the independent measure of the parameter pi, and p is the estimated value of the parameter pi. This augmented objective function is called the maximum a posteriori (MAP) estimator objective function. A recent application in physiological modeling is in Ref. 11.

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
Cov(<A><AC>p</AC><AC>ˆ</AC></A>)<IT>≥</IT><FENCE><FENCE><FR><NU><IT>∂</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t</IT>)</NU><DE><IT>∂</IT><A><AC>p</AC><AC>ˆ</AC></A></DE></FR></FENCE><SUP>T</SUP><IT>V<SUP>−1</SUP> </IT><FR><NU><IT>∂</IT>m(<A><AC>p</AC><AC>ˆ</AC></A><IT>, t</IT>)</NU><DE><IT>∂</IT><A><AC>p</AC><AC>ˆ</AC></A></DE></FR><IT>+&Sgr;<SUP>−1</SUP></IT></FENCE><SUP><IT>−1</IT></SUP> (10)
With 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 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
PCr<SUB><IT>0</IT></SUB><IT>=</IT><FR><NU><IT>V</IT><SUB>rev CK</SUB>ATP<SUB><IT>0</IT></SUB>TCr<IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB></NU><DE><IT>V</IT><SUB>rev CK</SUB>ATP<SUB><IT>0</IT></SUB><IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB><IT>+V</IT><SUB>for CK</SUB>ADP<SUB><IT>0</IT></SUB><IT>K</IT><SUB>iq</SUB><IT>K</IT><SUB>p</SUB></DE></FR> (11)
and
ADP<SUB><IT>0</IT></SUB><IT>=K</IT><SUB>ADP</SUB><FENCE><FR><NU><IT>k</IT><SUB>rest</SUB>ATP<SUB><IT>0</IT></SUB></NU><DE><IT>V</IT><SUB>max</SUB><IT>−k</IT>ATP<SUB><IT>0</IT></SUB></DE></FR></FENCE><SUP><IT>1/n</IT>H</SUP> (12)
Estimates 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)
V<SUB>rev CK</SUB><IT>=</IT><FR><NU><IT>V</IT><SUB>for CK</SUB></NU><DE><IT>K</IT><SUB>eq</SUB></DE></FR> <FENCE><FR><NU><IT>K</IT><SUB>iq</SUB><IT>K</IT><SUB>p</SUB></NU><DE><IT>K</IT><SUB>ia</SUB><IT>K</IT><SUB>b</SUB></DE></FR></FENCE> (13)
The initial pH is typically 7.0 to 7.1 in human muscle. The following relationship relates Keq of the creatine kinase reaction to pH measurements obtained during the course of the experiment
K<SUB>eq</SUB><IT>=1.77×10</IT><SUP>(<IT>9 − </IT>pH)</SUP> (14)

                              
View this table:
[in this window]
[in a new window]
 
Table 2.   Parameters and constants in Model 1

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 s-1 · 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.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 2.   Sample data set from Blei et al. (6), where phosphocreatine (PCr) decline during prolonged ischemia was monitored as a measure of basal metabolic rate. Circles denote measured PCr concentrations, whereas the smooth line is the fit of the basic model (Eq. 3).


                              
View this table:
[in this window]
[in a new window]
 
Table 3.   Estimates of basal metabolic rate (kbasal) from data originally presented in Fig. 3 of Blei et al. (6) for four subjects

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
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
REFERENCES

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.

                              
View this table:
[in this window]
[in a new window]
 
Table 4.   Results of Model 1 analysis

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


View larger version (37K):
[in this window]
[in a new window]
 
Fig. 3.   Average data and model fit for Model 1. Model prediction (continuous line), averaged over all subjects, is shown together with PCr and Pi average data (squares and circles, respectively, with SD error bars).

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.

                              
View this table:
[in this window]
[in a new window]
 
Table 5.   Results of Model 2 analysis

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.


View larger version (36K):
[in this window]
[in a new window]
 
Fig. 4.   Average data and model fit for Model 2. Model prediction (continuous line), averaged over all subjects, is shown together with PCr and Pi average data (squares and circles, respectively, with SD error bars).

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, phi , 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
OxPhos(<IT>t</IT>)<IT>=</IT><FR><NU><IT>&phgr;·V</IT><SUB>max</SUB><FENCE><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ADP</SUB></DE></FR></FENCE><SUP><IT>n</IT>H</SUP></NU><DE><IT>1+</IT><FENCE><FR><NU><IT>&Dgr;</IT>ADP(<IT>t</IT>)</NU><DE><IT>K</IT><SUB>ADP</SUB></DE></FR></FENCE><SUP><IT>n</IT>H</SUP></DE></FR> (15)
We chose phi  = 3 after preliminary simulations where phi  varied between 1 and 10. In the modified model, this feed-forward activation (increase in the value) of Vmax is added during the recovery phase (phase E) only, and not at rest (during phases A, B, and F). Parameter results are shown in Table 6. The quality of fit improved substantially (Fig. 5), primarily during the resting phases of the experiment (often leading to lower estimates of krest close to the independently measured value) and during the recovery phase. The model resulted in slightly lower estimates for KADP without major differences in kstim (contractile ATPase flux). Overall, the parameter values were similar to those found in the other two models. We also tested feed-forward increases in KADP with similar results; however, KADP was more sensitive to the magnitude of phi  than was Vmax, because the term for KADP has the exponent nH. 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.

                              
View this table:
[in this window]
[in a new window]
 
Table 6.   Results of Model 3 analysis



View larger version (35K):
[in this window]
[in a new window]
 
Fig. 5.   Average data and model fit for Model 3. Model prediction (continuous line), averaged over all subjects, is shown together with PCr and Pi average data (squares and circles, respectively, with SD error bars).


                              
View this table:
[in this window]
[in a new window]
 
Table 7.   AIC and BIC for the three models of cellular energy balance

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
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
REFERENCES

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.

                              
View this table:
[in this window]
[in a new window]
 
Table 8.   P values for the differences in AIC and BIC for Models 1, 2, and 3 

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
TOP
ABSTRACT
INTRODUCTION
EXPERIMENTAL DESIGN AND METHODS
RESULTS
DISCUSSION
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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

35.   Pette, D, and Staron RS. The molecular diversity of mammalian muscle fibers. News Physiol Sci 8: 153-157, 1993[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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[Abstract/Free Full Text].

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


Am J Physiol Cell Physiol 279(1):C213-C224
0363-6143/00 $5.00 Copyright © 2000 the American Physiological Society