Metabolic Dynamics in Skeletal Muscle during Acute Reduction in Blood Flow and Oxygen Supply to Mitochondria: In-Silico Studies Using a Multi-Scale, Top-Down Integrated Model

Control mechanisms of cellular metabolism and energetics in skeletal muscle that may become evident in response to physiological stresses such as reduction in blood flow and oxygen supply to mitochondria can be quantitatively understood using a multi-scale computational model. The analysis of dynamic responses from such a model can provide insights into mechanisms of metabolic regulation that may not be evident from experimental studies. For the purpose, a physiologically-based, multi-scale computational model of skeletal muscle cellular metabolism and energetics was developed to describe dynamic responses of key chemical species and reaction fluxes to muscle ischemia. The model, which incorporates key transport and metabolic processes and subcellular compartmentalization, is based on dynamic mass balances of 30 chemical species in both capillary blood and tissue cells (cytosol and mitochondria) domains. The reaction fluxes in cytosol and mitochondria are expressed in terms of a general phenomenological Michaelis-Menten equation involving the compartmentalized energy controller ratios ATP/ADP and NADH/NAD+. The unknown transport and reaction parameters in the model are estimated simultaneously by minimizing the differences between available in vivo experimental data on muscle ischemia and corresponding model outputs in coupled with the resting linear flux balance constraints using a robust, nonlinear, constrained-based, reduced gradient optimization algorithm. With the optimal parameter values, the model is able to simulate dynamic responses to reduced blood flow and oxygen supply to mitochondria associated with muscle ischemia of several key metabolite concentrations and metabolic fluxes in the subcellular cytosolic and mitochondrial compartments, some that can be measured and others that can not be measured with the current experimental techniques. The model can be applied to test complex hypotheses involving dynamic regulation of cellular metabolism and energetics in skeletal muscle during physiological stresses such as ischemia, hypoxia, and exercise.


Introduction
Skeletal muscle plays a major role in the regulation of wholebody substrates and energy metabolism, especially under changing physiological conditions such as ischemia (reduced blood flow), hypoxia (reduced oxygen supply), and exercise (increased energy demand). Current experimental techniques provide relatively little in vivo data on dynamic responses of muscle metabolite concentrations and metabolic fluxes to such physiological stimuli, especially in subcellular domains, such as mitochondria. To quantitatively analyze available in vivo experimental data and predict nonmeasurable dynamic responses, we developed a physiologically-based, multi-scale computational model of skeletal muscle cellular metabolism and energetics. The model is developed here from our previous model of cellular metabolism and energetics in skeletal muscle [1] and incorporates interdomain transport processes and compartmentalized metabolic reactions of many key chemical species in both cytosol and mitochondria.
Developing a mechanistic computational model of substrates and energy metabolism in complex, multi-scale metabolic systems, such as skeletal muscle, using a detailed, bottom-up systems approach with sparse in vivo experimental data-with an objective of achieving a quantitative understanding of the system to physiological perturbations-is a challenging task. Such a modeling approach requires information about the general structural features and catalytic mechanisms of the associated transporters and enzymes, subcellular metabolic pathways and fluxes and their control mechanisms, and tissue/organ specific metabolic characteristics. Such a modeling approach also requires mechanistic models for key functional components of the system (e.g., inter-domain transport processes, glycolysis, TCA cycle, oxidative phosphorylation, fatty acid b-oxidation) to be first individually developed and validated and then integrated to emulate the systems behavior at the molecular, subcellular, cellular, and tissue/organ levels. To avoid this complex approach and facilitate analysis of available sparse in vivo experimental data to understand dynamic responses of the system to physiological stresses, approximations are often made to obtain a simplified model of the system that includes key functional components regulating cellular metabolic processes at the desired level of complexity.
A top-down systems approach is an alternative approach which has been previously applied to determine and integrate a representative set of lumped biochemical reactions in vivo metabolic systems that incorporate primary substrates and key intermediate metabolites with coupled metabolic energy controllers ATP-ADP and NADH-NAD + [1][2][3][4][5][6][7][8]. This approach is similar to the top-down systems approach in metabolic control analysis proposed by Brand and co-workers [9,10] and is intended to provide an essential or minimal set of stoichiometrically balanced lumped biochemical reactions participating in ATP synthesis within mitochondria from the metabolism of nutrients (e.g., glucose, fatty acids, amino acids). Even with such simplifications, a large number of phenomenological kinetic parameters are introduced in the governing model equations, which must be estimated from available sparse in vivo experimental data. To estimate these unknown parameters, constraint-based, robust nonlinear optimization methods are needed, as established in our previous work [1].
To date, no physiologically-based, whole-organ level model of skeletal muscle cellular metabolism and energetics has been developed that can be applied to analyze available in vivo experimental data and predict dynamic metabolic responses to physiological stimuli in subcellular compartments, such as mitochondria. Previous models have incorporated some aspects of glycolysis, TCA cycle, oxidative phosphorylation, and fatty acid b-oxidation [3][4][5][11][12][13][14][15][16][17][18][19][20][21]. None of them, however, include sufficient key substrates/metabolites and/or integrate metabolic pathways/reactions that are essential in the regulation of cellular metabolic processes in skeletal muscle in vivo at whole-organ level. While our recently developed models of skeletal muscle cellular metabolism and energetics [1] have incorporated key metabolic pathways and reactions, the intracellular cytosolic and mitochondrial compartments were not distinguished. The energy controller pairs ATP-ADP and NADH-NAD + that modulate several key metabolic reactions in the cytosol and mitochondria have different concentrations in these two subcellular domains [6][7][8]. As a consequence, the model may not accurately predict the dynamics of several metabolite concentrations and metabolic fluxes that are critical in the regulation of fuel (carbohydrate, fat, and lactate) metabolism and cellular respiration during physiological stresses such as ischemia, hypoxia, and exercise.
In this paper, a physiologically-based, whole-organ level model of skeletal muscle cellular metabolism and energetics is developed and applied to study dynamic cellular metabolic responses to reduced blood flow and oxygen supply to mitochondria (muscle ischemia). The model, which is extended from our previous model [1], is based on a multi-scale, top-down systems approach [1][2][3][4][5][6][7][8], accounts for subcellular compartmentalization, and includes primary substrates (carbohydrates and fats) and key intermediate metabolites and metabolic reactions specific to skeletal muscle metabolic system. The model equations are based on dynamic mass balances of chemical species in capillary blood and tissue cells (cytosol and mitochondria) domains. The model also distinguishes the free and bound forms of O 2 and CO 2 transport in the blood and cells. The inter-domain species transport processes are considered either by passive diffusion or by carrier-mediated (facilitated) transport. The metabolic reaction fluxes in the cytosolic and mitochondrial domains are represented by a general phenomenological Michaelis-Menten equation involving the compartmentalized ATP/ADP and NADH/NAD + energy controller ratios. The phenomenological kinetic parameters of the model are estimated by using our recently developed constraintbased, robust nonlinear optimization approach [1]. In this estimation process, we fit the model output to sparse in vivo dynamic data on glycolytic and energetic metabolite concentrations from experiments in humans with muscle ischemia previously published [22]. With the estimated optimal parameter values, the model is able to simulate dynamic responses of key chemical species and reaction fluxes to reduced blood flow and oxygen supply to mitochondria associated with the muscle ischemia.

Model Development
The first step in our development of a multi-scale, top-down computational model of cellular metabolism and energetics in skeletal muscle was to identify the key intermediate metabolites and regulatory enzymes in the cellular metabolic pathways of skeletal muscle. We then integrated available information on cellular metabolic pathways and fluxes, cellular metabolic control mechanisms, catalytic enzyme kinetic mechanisms, subcellular compartmentation and metabolites volumes of distributions, interdomain transport mechanisms, and skeletal muscle tissue-specific metabolic characteristics. A simplified map of the compartmentalized cellular metabolic pathways in skeletal muscle is shown in Figure 1. The lumped biochemical reactions in the metabolic pathways are generated by stoichiometrically coupling several sequential elementary reactions. These reactions include the compartmentalized metabolic energy controller pairs ATP-ADP and NADH-NAD + whose ratios are known to modulate (fine tune) the reaction fluxes [2,23] in the subcellular compartments. Many of these lumped reactions are considered irreversible for which the resting Gibbs free energy (DG) is high and negative in favor of product formation [24]. As a part of a general formalism for modeling in vivo metabolic systems (nonequilibrium open systems), the reversible reactions like lactate dehydrogenase (LDH), creatine kinase (CK), and adenylate kinase (AK) were decomposed into two separate irreversible reactions with distinct kinetics [1].

Dynamic mass balance equations
The dynamic mass balance equations are based on a multidomain model structure for skeletal muscle consisting of a spatially-lumped capillary blood domain which exchanges nutrients and metabolic waste products with a spatially-lumped domain of tissue cells ( Figure 2). Although these two domains are separated by the interstitial fluid (ISF) space, we assume phase-equilibrium of chemical species between the blood and ISF domains, and consider them together as the ''blood'' domain. Furthermore, the The exchange arrows show the direction of net tissue cells uptake-release rates at normal, resting conditions. The lumped reactions are further compartmentalized into the cytosolic and mitochondrial reactions. These two subcellular domains are assumed to be in rapid equilibrium state (the barrier is shown schematically by double dotted lines), so that the species that are common to both these domains can have the similar dynamics, i.e., C mit,j (t) = s j .C cyt,j (t), where s j is the partition coefficient of the species j between cytosol and mitochondria. In this case, the net transport flux of a species j across the cytosol-mitochondria barrier can quickly become negligible at the onset of a physiological perturbation (the transport fluxes are shown by double dotted arrows); 10 species (PYR, FAC, CoA, NADH, NAD + , ATP, ADP, PI, CO 2 and O 2 ) exist in both the cytosolic and mitochondrial domains. GLC: glucose, GLY: glycogen, G6P: glucose-6-phosphate, GA3P: glyceraldehyde-3-phosphate, 13BPG: 1,3-biphosphate-glycerate, PYR: pyruvate, LAC: lactate, ALA: alanine, TGL: triglycerides, GLR: glycerol, FFA: free fatty acid, FAC: fatty acyl-CoA, ACoA: acetyl-CoA, CIT: citrate, AKG: a-ketogluterate, SCoA: succinyl-CoA, SUC: succinate, MAL: malate, OXA: oxaloacetate, CoA: coenzyme-A (free), PCR: phosphocreatine, CR: creatine, PI: inorganic phosphate, CO 2 : carbon dioxide, O 2 : oxygen, NADH: reduced tissue cells domain is compartmentalized into the cytosolic and mitochondrial domains. The chemical species are assumed to be distributed in these two subcellular domains as per their mass fractions and volumes of distributions (Table 1). Here, the mass fraction of a species in a particular domain is defined as the fractional amount of the species in that domain in comparison to the total amount of the species in the whole muscle tissue cells. The volume of distribution of a species in a particular domain is defined as the anatomical volume of the domain plus the binding space of the domain for the species. In addition, the species that are common to both of these domains are assumed to have the similar dynamics in these two domains, because the species transport processes between these two domains can be sufficiently fast [5]. Consequently, a change in the species concentration in one domain will be proportional to the change in the species concentration in the other domain: C mit,j (t) = s j .C cyt,j (t), where s j is the equilibrium concentration ratio (or partition coefficient) of species j between mitochondria and cytosol. In the present model, a total of 10 chemical species (pyruvate, fatty acyl-CoA, CoA, ATP, ADP, inorganic phosphate, NADH, NAD + , O 2 , and CO 2 ) are considered to exist in both the cytosolic and mitochondrial compartments with negligible transport flux between the compartments.
The dynamic mass balance of a chemical species j in the spatially-lumped blood domain has the following general form: where C art,j is the arterial species concentration; C bl,j is the capillary blood species concentration (equal to the venous species concentration C ven,j ); V bl,j and V isf,j are the volumes of distribution of species j in blood and ISF, and Q is the regional blood flow; J bl«cyt,j is the net transport flux (mass per unit time) across the blood-cytosol exchange barrier (consisting of capillary membrane, ISF, and tissue cell membrane). The dynamic mass balance of the chemical species j in the spatially-lumped tissue cells domain (cytosol and/or mitochondria) has the following general forms: where Here C x,j is the species concentration in domain x (cytosol/ mitochondria); V x,j is the volume of distribution of species j in domain x; P x,j and U x,j are the production and utilization of species . Schematic diagram of the structure of the model for blood-tissue cells exchange and cellular metabolism in skeletal muscle. The compartments are assumed to be perfectly mixed, and the capillary blood and tissue ISF regions are assumed to be in phaseequilibrium with each other, so that C isf,j = C bl,j = C ven,j for any chemical species j. The tissue cells domain is further compartmentalized into the cytosolic and mitochondrial domains with the chemical species having the similar dynamics in these two subcellular domains, so that C mit,j (t) = s j C cyt,j (t), where s j is the partition coefficient of the species j between cytosol and mitochondria. The model accounts for 30 chemical species in the tissue cells. A total of 8 species (GLC, LAC, PYR, ALA, FFA, GLR, CO 2 and O 2 ) undergo blood-tissue cells exchange; 10 species (PYR, FAC, CoA, NADH, NAD + , ATP, ADP, PI, CO 2 and O 2 ) exist in both the cytosolic and mitochondrial domains with a negligible transport flux (J cyt«mit,j <0). For details, see the caption of Figure 1 j in domain x; w x,p and w x,u are the reaction fluxes of the reactions processes that produce and utilize species j in domain x; b x,j,p and b x,j,u are the corresponding stoichiometric coefficients. For chemical species which are in the tissue cells (cytosol and/or mitochondria) but not in the capillary blood, the transport flux J bl«cyt,j is zero.
The dynamic mass balance equations for all the chemical species in blood, cytosol and mitochondrial domains can be rewritten from our previous model of skeletal muscle metabolism [1] by accounting for the species compartmentalized volumes of distributions as laid out in Eqs. (1) and (2a-2c). The dynamic mass balance equations for O 2 and CO 2 in these domains are developed by considering their distinct transport and binding mechanisms [25,26]. The detailed mass balance equations of the chemical species, including O 2 and CO 2 , are provided in Materials S1.

Transport and reaction flux equations
The reversible transport flux J bl«cyt,j (mass per unit time) of the species j across the blood-cytosol exchange barrier is related to the concentrations C bl,j and C cyt,j by The species concentrations (mM or mmol/kg tissue cells ww) are converted from the experimental data (mmol/kg tissue cells dw) in the literature by multiplying a factor 0.25 kg tissue cells dw/kg tissue cells ww [33]; for unit density, kg ww = L and mmol/kg ww = mmol/L = mM. The cytosolic and mitochondrial species concentrations are calculated from the species concentrations in the muscle tissue cells according to their approximate volumes of distributions and mass fractions; species nomenclature is adopted from Ref.
[1] (also see the caption to Figure 1); dw denotes dry weight, ww denotes wet weight. # The mass fractions are set to have reasonable compartmentalized species concentrations and mass action ratios consistent with available information from the literature: C cyt,LAC /C cyt,PYR = 16.4, C cyt,NAD+ /C cyt,NADH = 540, C mit,NAD+ /C mit,NADH = 6.3, C cyt,PCR /C cyt,CR = 2, C cyt,ATP /C cyt,ADP = 333.2, C mit,ATP /C mit,ADP = 1.11, K LDH = (C cyt,LAC / C cyt,PYR )* (C cyt,NAD+ /C cyt,NADH ) = 8856, K CK = (C cyt,CR /C cyt,PCR )*(C cyt,ATP /C cyt,ADP ) = 166.6, and K AK = (C cyt,ADP ) 2 /(C cyt,ATP C cyt,AMP ) = 1.2E-3. The cytosolic and mitochondrial PI concentrations are set at equal value with the assumption of a negligible pH gradient across the mitochondrial membrane. doi: 10 where l bl«cyt,j is the effective permeability surface area product for diffusive mass transport across the barrier (for passive transport); T bl«cyt,j is the maximal transport flux across the barrier (T max ) and M bl«cyt,j is the corresponding Michaelis-Menten (M-M) constant (M m ) (for facilitated transport). The flux expressions (3) satisfy the thermodynamic equilibrium conditions across the barrier. The species involved in the blood-cytosol exchange are glucose, lactate, pyruvate, alanine, glycerol, free fatty acid, CO 2 , and O 2 ( Table 2). The transport processes may be passive or carrier-mediated (facilitated) (Eq. 3). The transport of glucose, pyruvate/lactate, and free fatty acid across the sarcolemma of skeletal muscle is facilitated via the GLUT4 [27], MCT1 or MCT4 [28,29], and FABP or FAT/ CD36 [30] proteins, respectively. The sarcolemmal transport of the remaining species (alanine, glycerol, CO 2 , and O 2 ) is considered to be passive. For the 10 species that exist in both cytosol and mitochondria, the inter-domain transport flux J cyt«mit,j is zero.
The lumped metabolic reactions of skeletal muscle cellular metabolism and energetics can be considered as special cases of a general, irreversible, multi-reactant multi-product enzymatic reaction coupled with the metabolic energy controller pairs: where P 1 and P 2 are ATP and ADP or vice-versa (PS6: phosphorylation state); R 1 and R 2 are NADH and NAD + and vice-versa (RS6: redox state). The corresponding coupled phosphorylation and redox reactions are ATPRADP (PS+) and/or NADHRNAD + (RS+) or vice-versa (PS-and RS-). Assuming a phenomenological, single-step enzyme kinetic mechanism [31], the flux expression for the lumped metabolic reaction can be written as (see Ref. [1] for detailed description): where  . * A resting blood flow of 0.9 L/min corresponding to the two-leg quadriceps femoris muscle is used to calculate the muscle tissue cells uptake-release (UR) rates from the arterio-venous (AV) differences. The venous species concentrations and uptake-release rates are tuned further to satisfy the physiological constrains # stated above, just below the table. (These well-compiled data are adapted from ). Besides, there are 10 partition coefficients s for the chemical species that exists in both cytosol and mitochondria (since the 10 J cyt«mit,j transport fluxes are considered negligible). Therefore, the present model is characterized by a total of 99 unknown parameters that need to be estimated to fit the model outputs to the available in vivo experimental data. For convenience, the transport and reaction fluxes are written here in vector form:J J~J JC C;l l,T T max ,M M m and w w~w wC , where the parameter vector for these transport and reaction fluxes is denoted bỹ ;C C is the concentration vector. Note that the 12 transport parameters (4 l, 4 T max and 4 M m ) can be estimated here uniquely from the resting, steadystate transport flux-concentration relationships (Table 3) [1]. The 10 partition coefficients (10 s) can be estimated uniquely from the resting, steady-state species concentration ratios between the mitochondria and cytosol (s j = C mit,j /C cyt,j ) ( Table 3).

Model simulation
With specified parameter values, the mathematical model is solved numerically to simulate dynamic responses of the system to ischemia produced by reducing muscle blood flow. Typically, the initial conditions: C bl,j tƒ0 ð Þ~C 0 bl,j , C cyt,j tƒ0 ð Þ~C 0 cyt,j , and C mit,j tƒ0 ð Þ~C 0 mit,j are assumed to be at a normal, resting steady-state. These are fixed based on average resting species concentrations gathered from various literature sources on skeletal muscle cellular metabolism that are consistent with the resting, steady-state flux-concentration relationships (Tables 1 and 2). The species mass fractions and volumes of distributions in the subcellular cytosolic and mitochondrial domains are set to have appropriate phosphorylation and redox potentials in the cytosol and mitochondria. For numerical solution of this stiff initial-value problem, a robust implicit integrator DLSODES (https:// computation.llnl.gov/casc/odepack/odepack_home.html; http:// www.netlib.org/odepack; [32]) is used. Specifically, the Gear's implicit integration method based on backward difference formula (BDF) is most suitable for this problem. An absolute and relative error of tolerance of 10 210 guarantees high accuracy and convergence of the iterative solutions of the ODEs. The DLSODES solver is usually very fast; a typical simulation of this problem using the DLSODES solver in a standard desktop computer (Intel Xeon or Core 2 Duo CPU 5160 @ 3 GHz) takes only about 5 seconds of the CPU time. As a check, the numerical solutions of the initial value problem were also obtained using the ODE15S solver in MATLAB (http://www.mathworks.com) with Table 3. Optimal model parameter values for the inter-domain transport fluxes determined from the steady-state parameter estimation process.
I. Blood-cytosol passive or carrier-mediated (facilitated) transport fluxes (mmol/min) and parameters: Volumes of blood, ISF, cytosol and mitochondria Muscle blood flow at rest for two legs 0.9 L/min [33,57,59] Q isch Muscle blood flow during ischemia 0.216 L/min estimated The table also includes the anatomical volumes and muscle blood flow at rest (fixed) and during muscle ischemia (estimated).* * The blood-cytosol transport fluxes satisfies J bl«cyt,j = UR j = Q(C art,j 2C ven,j ) at resting, steady-state (the UR rates column in Table 2). The subcellular cytosolic and mitochondrial domains are assumed to be in rapid equilibrium state so that the cytosol-mitochondria transport fluxes are considered negligible and the metabolites common to both of these domains are considered to have the similar dynamics (i.e., C mit,j (t) = s j .C cyt,j (t)). Thus the cytosol-mitochondria metabolites partition coefficients s j are calculated based on the resting, steady-state metabolite concentrations in these two subcellular domains. The estimates of the transport flux parameters differ from those obtained in our previous work [1] because of the re-estimation of these parameters due to subcellular compartmentalization. doi:10.1371/journal.pone.0003168.t003 the similar tolerance levels. These solutions were of comparable accuracy with that obtained using the DLSODES solver.

Parameter Estimation
The large number of unknown parameters of this model are estimated by comparing model outputs to in vivo experimental data of Katz [22] with the optimization procedure described previously [1]. The data consist of key metabolites concentration dynamics measured during circulatory occlusion and recovery (reperfusion) in human skeletal muscle at the whole tissue-organ level [22]. Specifically, the data is based on the biopsy measurements of glucose and glycolytic/glycogenolytic intermediates (i.e., glucose 6phosphate, pyruvate, and lactate) and creatine and high-energy phosphates (i.e., phosphocreatine, inorganic phosphate, ATP, ADP, and AMP) in quadriceps femoris muscle tissue at rest, after 30 minutes of ischemia, and after 15 minutes of reperfusion. The tissue metabolite contents or concentrations were measured in the units of mmol/kg dry weight. For our analysis, these concentrations are converted to the units of mmol/kg wet weight (mmol/L or mM) by multiplying a conversion factor of 0.25 kg dry weight/ kg wet weight, corresponding to the muscle tissue [33]. Consistent with the experimental measurements of Katz [22], the following ischemia protocol is used for model simulations and parameter estimation: the muscle blood flow is reduced from Q = 0.9 L/min at rest (t,0 min) to Q = Q isch (unknown) at the onset of ischemia (0#t#30 min), and then increased to the starting level at the onset of recovery (t.30 min). The active muscle volume (V mus ) and ischemic muscle blood flow (Q isch ) responsible for the predicted metabolic dynamics in skeletal muscle during ischemia and recovery were not known from the experimental study of Katz [22]. Therefore, V mus and Q isch are also included as the unknown parameters for estimation, making the total number of unknown parameters in the model to 101. This is an increase of 10 unknown parameters (10 s) from our previous model [1].
With this relatively sparse data and large number of unknown parameters, the parameter estimation problem is ill-conditioned and under-determined. Nevertheless, an efficient estimation method is devised based on our earlier work [1] to obtain physiologically reasonable parameter values that minimize the sum of squared differences between the available experimental data and corresponding model outputs. The estimation procedure proceeds in two main stages. In the first stage, the published normal, resting species concentration data (Tables 1 and 2) are used to evaluate the resting transport and reaction fluxes from a steady-state flux balance analysis (Tables 2, 3 and 4). From the resting flux-concentration relationships: Eqs. (3) and (5), the preliminary estimates of the transport and reaction parameters are obtained. The preliminary estimates of the 10 partition coefficients are obtained from the resting, steady-state species concentration ratios between mitochondria and cytosol (s j = C mit,j /C cyt,j ) ( Table 3). In the second stage, the dynamic species concentration data [22] together with the resting, steady-state flux balance equations as equality constraints are used to obtain the optimal parameter estimates (Tables 3 and 4).
A detailed description of this robust estimation approach for obtaining optimal parameter estimates from species concentration dynamics during muscle ischemia and recovery is presented in Ref.
[1]. This method has been established powerful in the analysis of large-scale in vivo metabolic systems. The efficiency and robustness of this parameter estimation approach was tested with parameter sensitivity analysis as well as with repeated parameter estimation with various initial parameter estimates (initial guesses). Various estimates of the more sensitive model parameters spanned in a small neighborhood of the optimal parameter estimates (see Ref. [1] for details). Since the preliminary estimates of the 12 transport parameters and 10 partition coefficients based on the resting species concentrations and transport fluxes were accurate enough, these 22 parameters were not re-estimated from the dynamic species concentration data. Therefore, as in our previous paper [1], a total of 101-22 = 79 parameters (77 reaction parameters+V mus +Q isch ) are effectively estimated from the dynamic data.

Comparison of model simulations with experimental data
The optimal estimates of the transport parameters (l,T max ,M m ), partition coefficients (s), active muscle volume and ischemic muscle blood flow (V mus , Q isch ), and reaction parameters for the model are shown in Tables 3  and 4. These parameter estimates yield the best fit of the model outputs to the published experimental data [22] with minimal residual errors and minimal objective function. The optimal estimates of V mus and Q isch corresponding to the experimental data were ,4.0 L (,16% of normal two-legs muscle volume of ,25 L) and ,0.216 L/min (,76% reduction from normal, resting twolegs muscle blood flow of Q = 0.9 L/min). These optimal parameter values were used for model simulations during the resting, ischemia and recovery periods. The blood flow reduction levels of Q isch = 0.18, 0.27 and 0.36 L/min (80%, 70% and 60%) were used for simulating severe to moderate to mild ischemic conditions. The correspondence of model simulations and experimental data is demonstrated through Figures 3 and 4. Specifically, shown are the concentration dynamics of 4 glycolytic metabolites (GLC, G6P, LAC, PYR) and 6 energy metabolites (PCR, CR, PI, ATP, ADP, AMP) in the muscle tissue cells (i.e., weighted volume averages of concentrations in cytosol and mitochondria) for which experimental data were available [22] for four different blood flow reduction levels Q isch = 0.18, 0.216, 0.27 and 0.36 L/min.  Figure 5. The level Q isch = 0.216 L/min corresponds to the experimental data of Katz [22] and the corresponding model simulations match to the data reasonably well within the experimental noise ( Table 5). The experimental data are normalized here with respect to the resting metabolites concentrations (control). Furthermore, the individual metabolites responses are shown in separate plots in order to distinguish metabolite responses to different levels of blood flow reductions.
The model simulations of cellular glucose (GLC), glucose-6phosphate (G6P), pyruvate (PYR), and lactate (LAC) concentrations at the end of 76% ischemia and reperfusion are in close agreement with the experimental data ( Fig. 3(A-D), Table 5). Cellular [G6P] increased quickly (exponentially) by ,60% during ischemia and returned rapidly to its resting level at the onset of reperfusion. In contrast, the cellular [GLC] increased slowly (almost linearly) by ,35% during ischemia and remained at an elevated level even after 30 minutes of reperfusion. Blood [GLC], however, decreased rapidly, but only by ,10%, during ischemia and returned quickly to its baseline value during reperfusion (not shown). Furthermore, the model-predicted changes in the cellular [G6P] and [GLC] during the mild 60% and 70% blood flow reduction levels were not significant when compared to the changes during the high 76% and above (i.e., 80%) blood flow reduction levels.
Both ) decreased considerably from a resting level of ,8856 to ,5585 during 76% ischemia. At different levels of blood flow reductions, the model-predicted changes were not proportional. The changes were negligible below 70% blood flow reduction levels and significant above 76% blood flow reduction levels. The metabolic reaction flux rates are calculated based on resting, steady-state flux balance analysis (see Table 3, Ref. [1], for details). The flux rates with ''#'' can not be determined uniquely, and hence are set at values consistent with the other flux values (e.g., set a large value for the fast equilibrium reversible reactions like CK, AK and LDH). For K PS6 m and K RS6 m parameters, ''+'' indicates that the energy controller ratio is C ATP /C ADP or C ATP /C AMP or C NADH /C NAD+ , ''2'' indicates that the energy controller ratio is C ADP /C ATP or C AMP /C ATP or C NAD+ /C NADH , and ''0'' indicates that no controller ratio appears in the flux expression; ## denotes the estimated K m value K G6P GLC?G6P À Á for the inhibition of hexokinase reaction by G6P. doi:10.1371/journal.pone.0003168.t004 The model simulations of muscle phosphocreatine (PCR), creatine (CR), and inorganic phosphate (PI) concentrations at the end of 76% ischemia and recovery are in good agreement with the experimental data ( Fig. 4(A-C)). After 30 minutes, 76% ischemia resulted in ,17.5% decrease in muscle PCR content, which was fully resynthesized after 30 minutes of recovery. Muscle CR and PI contents increased by ,35% and ,130%, respectively, at the end of 76% ischemia and returned to their resting levels at the end of recovery.  (Fig. 4(D-F)). A 76% blood flow reduction resulted in

Simulated dynamics of mitochondrial metabolites
The model-simulated dynamic responses of mitochondrial FAC, ACoA, CIT, AKG, SCoA, SUC, MAL and OXA during the  ischemia and recovery periods with blood flow reduction levels of 80%, 76%, 70% and 60% (Q isch = 0.18, 0.216, 0.27 and 0.36 L/ min) are shown in Figure 6. Experimental data were not available for these key mitochondrial and TCA cycle intermediate metabolites. In fact, due to the limitations in the available experimental techniques, most of these subcellular metabolites can not be measured conveniently in skeletal muscle tissue cells in vivo at the whole tissue-organ level. However, the present compartmentalized model is able to simulate the dynamic responses of these subcellular metabolites to physiological stresses such as muscle ischemia. Model simulations show that the changes in these metabolites concentrations during muscle ischemia are not simply proportional to the extent of blood flow reduction. The changes are negligible during the mild 60% and moderate 70% blood flow reduction levels, but significant above the severe 76% blood flow reduction level.

Simulated dynamics of fluxes of exchangeable substrates
The model-simulated dynamic responses of transport and reaction rates associated with the blood-tissue cells exchangeable substrates (i.e., glucose, lactate, pyruvate, free fatty acid, CO 2 , and O 2 ) during the resting steady-state (t,0 min), ischemia (0#t#30 min), and recovery (t.30 min) periods with a blood flow reduction level of 76% (Q isch = 0.216 L/min) are shown in Figure 7. The corresponding capillary blood (venous blood) and tissue cells levels of O 2 and CO 2 are shown in Figure 8. The transport rates include the species net uptake/release rates (UR j = Q(C art,j 2C ven,j )) and net blood-tissue cells exchange rates J bl<cyt,j~J p bl<cyt,j or J f bl<cyt,j . The reaction rates include the species net production/utilization rates (R cl,j = P cl,j 2U cl,j ) in the tissue cells (cytosol and/or mitochondria). Most of these model predictions are similar to those from our previous model [1] which do not account for the intracellular compartmentalization. The magnitudes of changes with the present compartmentalized model, however, are lower and consistent with the lower estimate of the blood flow reduction level during muscle ischemia (present Q isch = 0.216 vs. previous Q isch = 0.135). Figure 7A shows that the net muscle glucose uptake (UR GLC ) changes sharply (decreases/increases) at the onset of ischemia/ recovery due to a step change (decrease/increase) in muscle blood flow, while the net glucose uptake by muscle tissue cells (J bl«cyt,GLC ) remains fairly constant during ischemia and recovery. This results in a fast adjustment (decrease/increase) in the venous glucose concentration and glucose A-V difference and a rapid recovery in muscle net glucose uptake. The net glucose utilization (U cl,GLC 2P cl,GLC ) through the hexokinase reaction decreases/ increases quickly during ischemia/recovery due to the increase/ decrease of [G6P] and G6P inhibition of hexokinase reaction. This is reflected in the predicted elevated cellular [GLC] during recovery (Fig. 3A). Several inhibition mechanisms of hexokinase reaction by G6P were tested. A general uncompetitive inhibition mechanism, which effectively modifies the K m as well as V max of the reaction [23,31], was found to be the most suitable for the model to fit to the data.
The net pyruvate UR and J bl«cyt rates are closely matched where both decrease/increase during ischemia/recovery, except at the onset of recovery, where UR PYR and J bl«cyt,PYR rates first sharply decrease and then exponentially increase to the baseline level (Fig. 7B, dashed and dotted lines). The pyruvate metabolism switches from a net utilization to a net production and a net uptake to a net release during ischemia and vice-versa during recovery (as depicted through the solid line in Fig. 7B). Similar to pyruvate, the net UR and J bl«cyt rates of lactate were approximately matched, through the 2J bl«cyt,LAC rate was slightly higher than the 2UR LAC rate during ischemia and vice-versa during recovery (Fig. 7C,  dashed and dotted lines). However, the net lactate release (2UR LAC and 2J bl«cyt,LAC ) first decreases when muscle blood flow is reduced and then increases steadily (linearly) slightly above the baseline value during the remaining ischemic period. During recovery, the opposite, but more pronounced, trend occurs. The net lactate production (P cl,LAC 2U cl,LAC ) increased quickly at the onset of ischemia and reached a new steady-state level within the first 10 minutes of ischemia. At the onset of recovery, the net lactate production decreased quickly to the baseline value (solid line, Fig. 7C). Figure 7D shows that the net UR and J bl«cyt rates of free fatty acid (FFA) decreased quickly at the onset of ischemia and then  Fig. 7(D,E)). However, the CO 2 production was closely correlated to the O 2 consumption (solid lines, Fig. 7(E,F)). The muscle O 2 uptake, cellular O 2 uptake, and cellular O 2 consumption were all reduced during the ischemia period and were closely matched (all reduced by ,22% during 76% ischemia), except at the onset of ischemia and recovery, where the muscle O 2 uptake decreased/increased sharply with the step change (decrease/increase) in muscle blood flow and returned to the new steady state levels as ischemia progressed (Fig. 7F). The O 2 partial pressure (P O2 ) in the muscle tissue cells decreased from 25 mmHg to ,0.5 mmHg (i.e., about the K m value for O 2 consumption), while the corresponding CO 2 partial pressure (P CO2 ) increased from 46 mmHg to ,52 mmHg during the 76% ischemia period (Fig. 8(B,D)). Accordingly, the total O 2 content in the capillary blood (venous blood) decreased from ,6.6 mM (,36 mmHg) to ,0.4 mM (,8 mmHg) and the total CO 2 content in the blood increased from ,25.5 mM (,43.6 mmHg) to ,29.6 mM (,51 mmHg) (Fig. 8(A,C)). These model simulations indicate that the muscle tissue cells were hypoxic during the 76% ischemia period; the reduced O 2 consumption was maintained by the reduced O 2 supply during ischemia.

Simulated dynamics of fluxes associated with ACoA and ATP
The model-simulated dynamic responses of metabolic fluxes associated with the production and utilization of ACoA and ATP during the resting steady-state (t,0 min), ischemia (0#t#30 min) and recovery (t.30 min) periods with a blood flow reduction level of 76% (Q isch = 0.216 L/min) are shown in Figure 9. Specifically, Figures 9A and 9B show the relative contributions from carbohydrates (PYR) and fats (FAC) to the production of ACoA in mitochondria, and the dynamics of ACoA utilization to CIT production in the mitochondrial TCA cycle. Figures 9C and 9D depict the utilization rate of ATP through ATP hydrolysis in cytosol and production rate of ATP through oxidative phosphorylation (O 2 consumption) in mitochondria and creatine kinase buffer reaction in cytosol. Figure 7. Model-predicted dynamic responses of transport and metabolic fluxes of glucose, pyruvate, lactate, free fatty acid, oxygen, and carbon dioxide during the resting, ischemia and recovery periods with a blood flow reduction of about 76%. The responses were computed using the estimated optimal parameter values with the ischemia protocol of 25 to 0 min of resting, 0 to 30 min of ischemia, and 30 to 60 min of recovery. The muscle blood flow Q is reduced as a step from 0.9 L/min at rest to Q isch = 0.216 L/min at the onset of ischemia and returned to 0.9 L/min at the onset of recovery. Q(C a -C v ): uptake-release rates, J b«c : blood-tissue cells transport rates, P c -U c : tissue cells production-utilization rates. doi:10.1371/journal.pone.0003168.g007 The production of mitochondrial ACoA from PYR decreased rapidly at the onset of ischemia and then increased relatively slowly to reach a new steady-state level as ischemia progressed (biphasic behavior). During recovery, the trend was opposite (Fig. 9A, dashed line). The production of mitochondrial ACoA from FAC varied relatively little during the entire ischemic and recovery periods (Fig. 9B, dotted line). An increase in total production of ACoA from PYR and FAC occurred at the beginning of ischemia due to a transient increase in the mitochondrial fatty acid b-oxidation. However, as ischemia progressed, ACoA production decreased and reached a new steady-state level significantly below its baseline value (Fig. 9B, solid line). In contrast, the utilization of ACoA to produce CIT decreased at the onset of ischemia and then increased to match the total ACoA production from PYR and FAC (Fig. 9B, dashed-dotted line). A mismatch between the total production and utilization of ACoA at the onset of ischemia is reflected in the predicted mitochondrial accumulation of ACoA during ischemia (Fig. 6B).
The rate of ATP utilization in cytosol through ATP hydrolysis slowly decreased (exponentially) by about 20% during 76% ischemia and returned to its baseline level during the recovery period (Fig. 9C, dashed line). However, the rate of ATP synthesis in mitochondria through oxidative phosphorylation (O 2 utilization) rapidly decreased (as a step) by about 22% during 76% ischemia (Fig. 9C, dotted line) due to a rapid reduction in O 2 delivery to mitochondria by capillary blood flow (Fig. 8B), in spite of sufficient increases in the mitochondrial redox potential [NADH]/[NAD + ] (Fig. 5B) and phosphorylation potential [ADP]/[ATP] (Fig. 5F) (the activators of the lumped oxidative phosphorylation reaction). This transient deficit in ATP (Fig. 9D solid line) was mostly matched by the ATP production from the creatine kinase ATP buffer reaction (Fig. 9D, dashed-dotted line). The ATP supply from glycolysis/glycogenolysis (ATP deficit -ATP supply from creatine kinase reaction) during ischemia and recovery was almost constant, indicating that glycolysis/glycogenolysis is not a major source of ATP supply during partial ischemia. The time to reach new steady states (response time) for the ATP hydrolysis and creatine kinase reactions were about 15 minutes, while that for the ATP synthesis (or O 2 consumption) reaction was about 5 minutes. Thus, in contrast to our previous model [1], the present compartmentalized model is able to correctly simulate the dynamic responses of cytosolic and mitochondrial metabolic reaction fluxes that are responsible in the production and utilization of ACoA and ATP. Furthermore, the model is able to link substrate metabolism (carbohydrates vs. fat) to energy metabolism in skeletal muscle during physiological perturbations such as muscle ischemia.

Discussion
Modeling framework for analysis of sparse in vivo experimental data Analysis of dynamic changes in cellular metabolism and energetics in skeletal muscle in vivo in response to reduced blood   Figure 8. Model-predicted dynamic responses of O 2 and CO 2 (total contents and partial pressures) in the capillary blood (or venous blood) and tissue cells during the resting, ischemia and recovery periods with a blood flow reduction of about 76%. The responses were computed using the estimated optimal parameter values with the ischemia protocol: if (t,0 min | t.30 min) Q = 0.9 L/min, else Q = Q isch = 0.216 L/min. The total O 2 and CO 2 contents are based on various forms of O 2 and CO 2 transports in the capillary blood and tissue cells (see Materials S1). doi:10.1371/journal.pone.0003168.g008 flow and oxygen supply to mitochondria (ischemia) is limited by available in vivo experimental data. Relatively few in vivo measurements can be obtained of metabolite concentrations and metabolic fluxes in subcellular compartments, such as mitochondria. As a complementary approach to experimental measurements and as a framework for quantitatively analyzing available in vivo data, a physiologically-based, whole-organ level model of skeletal muscle cellular metabolism and energetics is developed. This model emphasizes a multi-scale, top-down systems approach [1-8] for quantitative understanding of metabolic processes at the molecular, subcellular and cellular levels in relation to the tissueorgan responses. For this purpose, integration of information was required from cellular metabolic pathways and fluxes of substrate and energy metabolism in cytosol and mitochondria, cellular metabolic control mechanisms, catalytic enzyme kinetic mechanisms, subcellular compartmentation and metabolites volumes of distribution, inter-domain transport mechanisms, and tissuespecific skeletal muscle metabolic characteristics. Since the extent to which enzyme activities and kinetic (Michaelis-Menten) constants determined from in vitro experimental studies are applicable to in vivo conditions is uncertain, our multi-scale, topdown modeling approach provides a unified, systematic framework for modeling and analysis of in vivo complex metabolic systems. A major improvement of the present model over our  (Fig. 1). The governing model equations are based on dynamic mass balances of chemical species in spatially-lumped, capillary bloodinterstitial fluid domain and two distinct intracellular cytosolic and mitochondrial domains (Fig. 2). These include the dynamic mass balances of O 2 and CO 2 which are developed by accounting for their distinct transport and binding mechanisms (hemoglobin and myoglobin-mediated) in these three domains [25,26]. However, since not much information is available regarding the blood flow  heterogeneity and fiber type distributions in skeletal muscle, these metabolic characteristics were not considered in the present study. Instead, different fiber types were lumped together to an effective active muscle volume with an effective blood flow supplying oxygen and other nutrients to the lumped muscle volume, which is sufficient to investigate the dynamic metabolic responses to muscle ischemia (also see discussions in Section 4.4). Furthermore, since the phenomenon of muscle acidification/alkalization (i.e., proton handling) under physiological stresses, such as muscle ischemia, is highly complex [17,21], this was not accounted for in the present model.
In this model, the compartmentalized lumped biochemical reactions along the cellular metabolic pathways are the stoichiometrically coupled sequential elementary reactions, which include the compartmentalized metabolic energy controller pairs ATP-ADP and NADH-NAD + whose ratios (phosphorylation and redox potentials) modulate (or fine tune) the reaction fluxes in the subcellular domains [2,23]. The lumped reactions for which the resting Gibbs free energy (DG) is large and negative in favor of product formation are considered irreversible [24]. The reversible reactions like lactate dehydrogenase (LDH), creatine kinase (CK), and adenylate kinase (AK) are decomposed into two separate irreversible reactions with distinct kinetics. This can be justifiable based on the fact that the in vivo metabolic systems are nonequilibrium open systems and the biochemical reactions usually operate far from true equilibrium under physiological stimuli (also see discussions in Section 4.3). Besides, the consideration of the law of microscopic reversibility (thermodynamics) requires a detailed knowledge of proton handling in cytosol and mitochondria [17,21], which is not considered in this phenomenological modeling study.  [23,31]. In comparison to the previous models in the literature [3][4][5][11][12][13][14][15][16][17][18][19][20][21], the present model is a selfconsistent, physiologically-based, multi-scale computational model dealing with cellular metabolism and energetics in skeletal muscle in vivo at the whole-organ level.
Optimal parameter estimation for the large-scale metabolic model A major challenge of this work was to evaluate a large number of unknown model parameters (precisely a total of 101) with relatively sparse experimental data from in vivo studies on muscle ischemia [22]. To deal with this under-determined, ill-conditioned problem, we had to devise a special optimization strategy (see Ref. [1] for details). Briefly, preliminary estimates of the model parameters were obtained based on resting, steady-state flux balance analysis together with the information on resting, steadystate metabolites blood and tissue cells (cytosol and mitochondria) concentrations, metabolites uptake-release (UR) rates or arteriovenous (AV) differences, and metabolic reaction flux rates gathered from various literature sources on skeletal muscle cellular metabolism in normal human subjects (Tables 1 and 2). With these preliminary parameter estimates, the model could not simulate the dynamic cellular metabolic responses to muscle ischemia. However, these estimates provided a good starting point in an iterative process of optimization to obtain the optimal parameter estimates (Tables 3 and 4) that yield the best least-squares fit of the model to the data [22] (Figs. 3 and 4). A good starting point is essential for optimization of large-scale problems with gradienttypes of optimization algorithms, such as the generalized reduced gradient algorithm (GRG2) [34], used in the present study. The first step of obtaining the preliminary parameter estimates provided accurate optimal estimates of the 12 transport parameters (4 l, 4 T max and 4 M m ) and 10 partition coefficients (10 s), resulting in a reduction of 22 kinetic parameters for estimation from the dynamic experimental data from muscle ischemia.
The key to success of our parameter estimation approach for obtaining the remaining 79 kinetic parameters (second step) was the incorporation of constraint information to the maximum extent in the GRG2 optimization algorithm [34]. These include (1) the bound constraints on kinetic parameters from experimental studies, (2) the nonlinear equality (physiological) constraints relating fluxes and concentrations at resting, steady-state, and (3) the nonlinear inequality (thermodynamic) constraints relating fluxes and Gibbs free energy for reversible reactions: w net .DG#0 [35][36][37]. Also, the bounded estimation of active muscle volume (V mus ) and ischemic blood flow (Q isch ) (not known from the experimental protocol of Katz [22]) provided additional degrees of freedom for fitting the model to the data. Use of such constraints reduced the number of unknown model parameters for simultaneous estimation (e.g., the 26 V max parameters could be explicitly estimated from the flux-concentration relationships at resting, steady-state). Such constraints also reduced uncertainty and enhanced accuracy in the estimated parameter values.
Since the available experimental data [22] were relatively sparse in comparison to the number of unknown model parameters, the parameter estimates for this under-determined, ill-conditioned problem were expected to be non-unique. The efficiency and robustness of this parameter estimation approach was tested previously [1] in detailed with repeated parameter estimations with various initial parameter estimates (guesses) and with parameter sensitivity analysis. Various parameter estimates of the more sensitivity model parameters spanned in a small neighborhood of the accepted optimal parameter estimates. Reestimation of the model parameters for the present extended compartmentalized model with the same parameter estimation approach and the quality of the model fitting to the experimental data (Figs 3 and 4) further testifies the robustness of the parameter estimation approach.

Model predictions related to muscle ischemia
This physiologically-based, multi-scale computational model of skeletal muscle cellular metabolism and energetics can simulate dynamic metabolic changes in cellular and subcellular compartments and provide quantitative analysis of the mechanisms of metabolic regulation during physiological stresses (e.g., reduced blood flow and oxygen supply to mitochondria associated with muscle ischemia). Model simulations based on the estimated optimal parameter values correspond well to the experimental data on muscle ischemia of Katz [22] (Figs. 3 and 4). The optimal estimates of active muscle volume V mus <4 L and ischemic muscle blood flow Q isch <0.216 L/min (Table 3) indicate that only ,16% of whole muscle volume (,25 L) was actively participating in cellular metabolism during reduced oxygen delivery to mitochondria induced by a blood flow reduction of ,76% from the resting level of 0.9 L/min. Though these estimates are slightly different from those estimated previously [1] (V mus <5 L and Q isch <0.135 L/min), when scaled to one leg with a resting blood flow of 0.45 L/min, the estimated V mus value becomes ,2 L which lies well within the experimental range of 1.5-3 L for one leg quadriceps femoris muscle [38][39][40][41]. The experimental data [22] used in this computational study were also based on measurements from similar muscle types. This indicates that the whole muscle is not uniformly perfused and that only a fraction of the muscle mass, in accord with metabolic demand, participates in cellular metabolic processes.
The  (Fig. 4). One of the key aspects in achieving these successful model simulations and good fittings was to model each of the reversible reactions like creatine kinase (CK) and adenylate kinase (AK) as two separate irreversible reactions with distinct kinetics. The forward and reversible reaction fluxes were considered regulated by the phosphorylation potentials [ATP]/[ADP] or [ADP]/[ATP] in the cytosol (Materials S2). Without considering such control mechanisms, it was not possible to simulate the experimental data on these high-energy metabolites during muscle ischemia and reperfusion.
We tested fitting the model to the data using alternative flux expressions for the CK and AK (also LDH) reversible reactions satisfying the Haldane relationship and thermodynamic equilibrium conditions. However, we were unable to achieve good fittings for the high-energy metabolites. This may be due to the limitations in the present model, or may be the fact that the data do not support thermodynamic equilibrium of the CK reaction during ischemia/recovery. It is evident from the data of Katz [22]  This issue may be resolved by including the dynamics and buffering of cytosolic and mitochondrial pH [17,21] and, to some extent, by distinguishing the dynamics of common metabolites within the subcellular cytosolic and mitochondrial compartments. The later will necessitate incorporation of transport mechanisms of common metabolites across the mitochondrial inner membrane, e.g., proton pumps, ATP/ADP translocase, PI/H + cotransporter, and exchange of TCA cycle substrates/intermediates [18,21]. Since ,98% of ADP is in mitochondria and only ,2% of ADP is in cytosol at normal, resting steady-state conditions, a small This information will be the mechanistic basis for the extension of the current model which can be used to analyze the data on highenergy phosphate metabolites.
With the estimated optimal parameter values, the present model is able to correctly simulate the dynamic metabolic responses of compartmentalized cytosolic and mitochondrial phosphorylation and redox potentials (Fig. 5), key mitochondrial and TCA cycle metabolite concentrations (Fig. 6), and other key metabolite concentrations and metabolic fluxes (Figs. 7-9) in skeletal muscle during ischemia and reperfusion for which no experimental data are available. Some of these simulations may be questionable without further validation, such as the predictions of cytosolic and mitochondrial phosphorylation and redox potentials (Fig. 5) Our compartmentalized model simulations show that the metabolic changes were negligible with respect to the extent of blood flow reduction up to ,70%, beyond which a profound derangement in substrate and energy metabolism occurred due to reduced supply of oxygen to mitochondria. The metabolic changes during severe ischemia (,76% and beyond) include (1) a switch from a net pyruvate uptake and utilization to a net pyruvate production and release (Figs. 3C and 7B); (2) an increase in the redox potential [NADH]/[NAD + ] (Figs. 3(E,F) and 5(A,B)) accompanied by a net formation and accumulation of cytosolic LAC and PYR ( Fig. 3(C,D)) and mitochondrial ACoA (Figs. 6B and 9(A,B)); (3) an increase in the glycolytic and creatine kinase contributions to ATP formation ( Fig. 9(C,D)) accompanied by a net depletion of cytosolic PCR (ATP stores) (Fig. 4A). In general, these predictions are in agreement with experimental results from studies under similar conditions [42][43][44]. However, since our model does not differentiate between FADH 2 and NADH, the less efficient fuel (i.e., fats) -from the oxygen consumption point of view-is considered to have the same oxygen cost as carbohydrates for ATP production.
However, this assumption may not be valid for the [ATP]/[ADP] and [NADH]/[NAD + ] ratios during ischemia/recovery, which are major modulators of several key metabolic reactions in cytosol and mitochondria. Therefore, the model may not accurately predict the dynamics of metabolite concentrations and metabolic fluxes that are critical in the regulation of fuel (carbohydrate, fat, and lactate) metabolism and cellular respiration during physiological perturbations, such as ischemia, hypoxia, and exercise. Furthermore, the inclusion of FADH 2 as a reducing equivalent (produced through oxidation of fatty acids) is essential when investigating the relative contribution of carbohydrate and fat oxidations in the synthesis of ATP in mitochondria.
Another limitation in the present model is that the model does not account for blood flow heterogeneity and fiber type distributions associated with skeletal muscle cellular metabolism. Specifically, the capillary blood domain is considered spatially-lumped and the advective-diffusive transport of chemical species in capillary blood is ignored. Such compartmental modeling may not account for situations where important gradients of chemical species exist (e.g., due to distinct metabolic characteristics associated with different fiber types). In skeletal muscle, where the blood flow is highly heterogeneous due to complex capillarization and fiber type distributions, more than 70% of arterial oxygen is extracted from the blood before it leaves the microcirculation [45]. High extraction results in large intracellular and intravascular gradients of oxygen and other nutrients. Therefore, the well-mixed assumption may not be valid and compartmental modeling may not account for the heterogeneity of oxygen transport and consumption in different fiber types. These issues can be addressed through spatially-distributed modeling of oxygen transport and consumption, coupled with other substrates and CO 2 [19,26].
Finally, in the present model most of the lumped biochemical reactions in the cellular metabolic pathways are considered irreversible in the direction of product formation. Although the flux in the reverse direction is almost one order of magnitude smaller than the flux in the forward direction under normal, resting steady-state conditions, their magnitudes may change during physiological perturbations such as ischemia, hypoxia or exercise. Therefore, it would be advantageous to consider all the biochemical reactions to be reversible and their fluxes to satisfy the thermodynamic equilibrium constraints [37], including the dynamics and buffering of pH [17,21], which can affect the apparent equilibrium constants and standard Gibbs free energy of the reactions.

Future model developments and potential applications
Including the substrate and cation transport mechanisms between cytosol and mitochondria and distinguishing common metabolites dynamics within these subcellular domains would be important in the analysis of dynamic cellular metabolic responses to high intensity exercise. The subcellular compartmentalization could account for distinct volumes of metabolites distribution, for example, metabolic channeling for glycolysis [6][7][8]. This is a way to recognize that the key glycolytic enzymes are bound together to form a multi-enzyme complex near the sarcolemma and sarcoplasmic reticulum. To investigate the effects of spatial locality due to blood flow heterogeneity and fiber type distributions on metabolic responses to muscle ischemia, hypoxia or exercise, one will require spatially-distributed models involving parallally-arranged capillaries representing different muscle fiber types. Furthermore, the use of alternative kinetic flux expressions based on a Michaelis-Menten formalism for reversible enzymatic reactions would satisfy the thermodynamic equilibrium conditions and would provide additional thermodynamic constraints for kinetic parameter estimation [37]. As in our previous paper [1], a formal dynamic parameter sensitivity analysis can be carried out to understand the relative importance of the model parameters on the measured outputs. This information can be used to improve confidence in the estimates of key model parameters by fixing the values of the least sensitive parameters (parameter space reduction approach).
With these enhancements and with additional metabolic pathways similar to those of the models of cellular metabolism and energetics in cardiac muscle [6][7][8], such as the inclusion of FADH 2 as a reducing equivalent, the advanced model of cellular metabolism and energetics in skeletal muscle can be used to analyze physiological responses to exercise (increased energy demand). Specifically, the model can assist in providing insight into long-standing or paradoxical questions, such as the extent of control that mitochondrial oxygen concentration and redox state exert over the rate of lactate production in the cytosol during heavy intensity exercise. In general, a multi-scale computational model of cellular metabolism and energetics in tissue/organ systems within the body that incorporate sufficient mechanistic processes can be used to (1) analyze and interpret in vivo experimental data; (2) provide new hypotheses as well as suggest how to test a given hypothesis; (3) design critical experiments; (4) quantify and predict dynamic responses to physiological stimuli that can not be directly measured; (5) evaluate relative importance of metabolic pathways and fluxes and their regulatory mechanisms under both normal and pathological conditions; and (6) provide the mechanistic basis for simulating integrated effects of altering enzyme activities and substrate concentrations through pharmacological agents and/or dietary inputs.

Supporting Information
Materials S1