Mechanistic analysis of multi-omics datasets to generate kinetic parameters for constraint-based metabolic models
- Cameron Cotten^{1, 2} and
- Jennifer L Reed^{1, 2}Email author
https://doi.org/10.1186/1471-2105-14-32
© Cotten and Reed; licensee BioMed Central Ltd. 2013
Received: 18 April 2012
Accepted: 16 January 2013
Published: 30 January 2013
Abstract
Background
Constraint-based modeling uses mass balances, flux capacity, and reaction directionality constraints to predict fluxes through metabolism. Although transcriptional regulation and thermodynamic constraints have been integrated into constraint-based modeling, kinetic rate laws have not been extensively used.
Results
In this study, an in vivo kinetic parameter estimation problem was formulated and solved using multi-omic data sets for Escherichia coli. To narrow the confidence intervals for kinetic parameters, a series of kinetic model simplifications were made, resulting in fewer kinetic parameters than the full kinetic model. These new parameter values are able to account for flux and concentration data from 20 different experimental conditions used in our training dataset. Concentration estimates from the simplified kinetic model were within one standard deviation for 92.7% of the 790 experimental measurements in the training set. Gibbs free energy changes of reaction were calculated to identify reactions that were often operating close to or far from equilibrium. In addition, enzymes whose activities were positively or negatively influenced by metabolite concentrations were also identified. The kinetic model was then used to calculate the maximum and minimum possible flux values for individual reactions from independent metabolite and enzyme concentration data that were not used to estimate parameter values. Incorporating these kinetically-derived flux limits into the constraint-based metabolic model improved predictions for uptake and secretion rates and intracellular fluxes in constraint-based models of central metabolism.
Conclusions
This study has produced a method for in vivo kinetic parameter estimation and identified strategies and outcomes of kinetic model simplification. We also have illustrated how kinetic constraints can be used to improve constraint-based model predictions for intracellular fluxes and biomass yield and identify potential metabolic limitations through the integrated analysis of multi-omics datasets.
Keywords
Background
Constraint-based models consider the physico-chemical constraints that exist on metabolism, including mass balances, flux capacities, thermodynamics, and transcriptional regulation of metabolic enzymes [1, 2]. One common constraint-based modeling approach, flux balance analysis (FBA) uses mass balance and reaction reversibility constraints to predict metabolic fluxes that maximize flux through a reaction or combination of reactions [2]. FBA typically is used to find flux distributions that maximize biomass production or ATP generation per total flux [3]. Some constraint-based models add transcriptional regulatory constraints to further limit flux values [4-6] when the associated enzymes are known to be unexpressed in certain conditions. Thermodynamic constraints have also been imposed to limit the direction a given reaction can proceed, and thermodynamic metabolic flux analysis (TMFA) uses these constraints to ensure that reactions only proceed in the thermodynamically-favorable direction [7-10]. TMFA models were some of the first constraint-based models that directly accounted for intracellular metabolite concentrations.
Constraint-based models and kinetic differential equation models have largely been divergent methodologies in systems biology. In constraint-based modeling, omission of kinetic considerations is generally seen as an advantage of the methodology, since determining kinetic information for an entire metabolic network is currently infeasible. However, FBA and TMFA predictions are not always consistent with experimental observations because of kinetic limitations on native enzymes. To account for these kinetic limitations, there have been some efforts to integrate kinetic expressions into constraint-based models. Beg et al.[11] were able to improve predictions of cellular growth rates by constraining fluxes using individual enzyme activity, enzyme volume needed to achieve a given flux, and the total enzyme volume. More recently, Yizhak et al.[12] integrated metabolomic and proteomic data into Michaelis-Menten kinetic expressions using in vitro K_{ m } parameters to more accurately predict flux changes (increases/decreases) than methods that did not consider kinetics. Similarly, dynamic flux balance analysis models have also used Michaelis-Menten kinetic expressions to constrain substrate uptake rates based on reactor concentrations [13-15].
In the present study, we estimated kinetic parameters in a kinetic model of central Escherichia coli metabolism by integrating fluxomic, proteomic, and metabolomics data. Data published by Ishii et al.[16] from single deletion mutants and the parental strain (at different growth rates) were used to construct a weighted sum of least squares (WSLS) parameter estimation problem using kinetic rate laws. Initially rate laws reported by Chassagnole et al.[17] were used; however, these rate laws resulted in parameters with large confidence intervals. A simplified kinetic model was subsequently constructed that resulted in smaller confidence intervals for kinetic parameters. Predictions from the kinetic model for metabolite concentrations and kinetic parameters were used to draw conclusions about metabolite-mediated inhibition and activation effects on enzymes, and distance from equilibrium for reactions in central metabolism. Using independent data sets, we then used the kinetic model to predict flux ranges that are consistent with estimated kinetic parameters and concentration data. We subsequently incorporated these flux ranges as constraints into a constraint-based model to improve predictions over FBA. This new constraint-based model with kinetic constraints is one of the most detailed constraint-based models with kinetic constraints to date, and is one of a few to only use in vivo estimates of kinetic parameters.
Methods
Constraint-based metabolic model
A central E. coli constraint-based metabolic model previously reported by Palsson [18] was used that included glycolysis, pentose phosphate, oxidative phosphorylation, and the citric acid pathways. This model was chosen because it includes the fluxes, proteins, and metabolites measured by Ishii et al.[16]. This model included mass balance constraints for central metabolic intermediates, as well as energy and redox carriers. Fluxes from Ishii et al.[16] were adjusted slightly to satisfy mass balance constraints and energy requirements in the metabolic models by minimizing the Euclidean distance between the fluxes in the constraint-based model and the reported flux distributions. These adjusted flux measurements were then used to estimate kinetic parameters and evaluate model-predicted fluxes from a constraint-based model with kinetic constraints. The most recent genome-scale model of E. coli (iJO1366) [19] was also used to evaluate the effects of kinetic constraints. A few reactions were excluded from the model (EDD, HEX1, F6PA, FBA3, FLDR2, and DRPA) so that the adjusted fluxes more closely matched the measured flux patterns.
Kinetic parameter estimation using multi-omic data
Here, I is the set of metabolites, K is the set of reactions with kinetic rate laws, L is the set of experimental conditions, and θ is a vector of all kinetic parameter values. Parameters except for equilibrium constants were assigned global upper and lower bounds of 10^{4} and 10^{-6}, respectively (in their respective units). The estimates of metabolite and enzyme concentrations (c_{ il } and e_{ kl }) are variables and are adjusted to be as close as possible to the experimentally-measured metabolite and enzyme concentrations (c_{ il }^{exp} and e_{ kl }^{exp}). Fluxes with kinetic rate laws (v_{ kl }) were fixed to the adjusted values from the constraint-based model. Biologically-relevant global upper and lower bounds were enforced on metabolite concentrations (0.001 mM and 10 mM) and enzyme concentrations (0.001 mg protein/gDCW and 10 mg protein/gDCW). Kinetic parameters were allowed to vary freely, except for equilibrium constants that were required to be within 35% to 285% of their in vitro measured values in standard conditions [21], corresponding to changes of 2.6 kJ/mol in the Gibbs free energy change of reaction. The weights in the WSLS objective function, W_{ il }^{ c } and W_{ kl }^{ e }, were assigned as the inverse of experimental variances for the metabolic and enzyme concentration measurements, respectively. This weighting method was chosen to keep residuals on the same scale and assign less importance to concentrations with measurements that were more uncertain. Experimental sample variances were calculated from concentration measurements from 4 biological replicates of the parental strain growing steadily at 0.2 h^{-1}.
The WSLS parameter estimation problem was solved using CONOPT3 as called by GAMS 23.3 (GAMS Development Corporation, Washington, DC) from multiple starting points to explore the non-convex feasible space and find multiple local optimal solutions. The feasible solution with the lowest objective value was selected as the set of optimal kinetic parameters values and concentrations. Confidence intervals for the set of kinetic parameters were estimated by determining the sensitivity of enzyme and metabolite concentrations to parameter values. This sensitivity was determined by making small perturbations to parameter values one at a time and re-optimizing Equations 1a-1e. These sensitivities were used to determine confidence intervals in a method described by Antoniewicz et al.[22]. Generally, smaller WSLS values and higher sensitivities lead to smaller confidence intervals and more confidence in parameters. Five conditions (Δpgm, ΔgpmA, Δzwf, ΔtktB, and parental strain at D=0.5 h^{-1}) were randomly drawn from the set of experimental conditions and left out during parameter estimation, to form an independent test set to later evaluate the resulting constraint-based model.
Predictions using FBA with kinetic constraints (KFBA)
Here, J is the set of all reactions, K is the subset of reactions with kinetic rate laws, and the binary variables y^{ + } and y^{ - } indicate whether an upper or lower kinetic bound is violated. A general lower bound (LB=0 or −1000 mmol/gDW/h for irreversible and reversible reactions, respectively) and an upper bound (UB=1000 mmol/gDW/h) were used for all reactions, including those with kinetic rate laws. All fermentative pathways were blocked (i.e., LB and UB set to 0), as only carbon dioxide was produced in experiments. The choice for UB and LB did not affect predictions, as fluxes did not go to these limits at the optimal solutions and the kinetic bounds were between the LB and UB values. The cellular growth rates (v_{biomass}) were fixed to the chemostat dilution rate. An analogous FBA problem was formulated using Equations 1a-2c. Since the FBA and KFBA solutions are not necessarily unique, we selected the alternate optimal solution with the lowest sum of squared flux values [23] by solving an additional minimization problem, which gives a unique answer. Residuals were calculated for the subset of measured fluxes as the squared difference between the adjusted flux measurement and FBA or KFBA flux prediction. When calculating the mean residual, only one flux residual was used for fluxes that must be directly proportional to one another based on the network stoichiometry (e.g., reactions in a linear pathway).
Results
In vivo kinetic parameter estimation
Rate laws used in the simplified model
v _{ PTS } | ${k}_{\mathit{cat}}{e}_{\mathit{PTS}}\frac{{c}_{\mathit{pep}}}{{c}_{\mathit{pyr}}}$ | v _{ PGI } | ${k}_{\mathit{cat}}{e}_{\mathit{PGI}}\left(1-\frac{{c}_{f6p}}{{c}_{g6p}{K}_{\mathit{eq}}}\right)$ |
v _{ PFK } | k _{ cat } e _{ PFK } C _{ atp } | v _{ ALDO } | $\frac{{k}_{\mathit{cat}}{e}_{\mathit{ALDO}}\left({c}_{\mathit{fdp}}-\frac{{c}_{\mathit{gap}}{c}_{\mathit{dhap}}}{{K}_{\mathit{eq}}}\right)}{{K}_{\mathit{fdp}}+{C}_{\mathit{fdp}}}$ |
v _{ TPI } | ${k}_{\mathit{cat}}{e}_{\mathit{TPI}}\left(1-\frac{{c}_{\mathit{gap}}}{{c}_{\mathit{dhap}}{K}_{\mathit{eq}}}\right)$ | v _{ GAPD } | ${k}_{\mathit{cat}}{e}_{\mathit{GAPD}}\left({c}_{\mathit{nad}}{c}_{\mathit{gap}}-\frac{{c}_{13\mathit{dpg}}{c}_{\mathit{nadh}}}{{K}_{\mathit{eq}}}\right)$ |
v _{ PGK } | ${k}_{\mathit{cat}}{e}_{\mathit{PGK}}\left(1-\frac{{c}_{\mathit{atp}}{c}_{3\mathit{pg}}}{{c}_{\mathit{adp}}{c}_{13\mathit{dpg}}{K}_{\mathit{eq}}}\right)$ | v _{ PGM } | ${k}_{\mathit{cat}}{e}_{\mathit{PGM}}\left(1-\frac{{c}_{2\mathit{pg}}}{{c}_{3\mathit{pg}}{K}_{\mathit{eq}}}\right)$ |
v _{ ENO } | ${k}_{\mathit{cat}}{e}_{\mathit{ENO}}\left({c}_{2\mathit{pg}}-\frac{{c}_{\mathit{pep}}}{{K}_{\mathit{eq}}}\right)$ | v _{ PYK } | ${k}_{\mathit{cat}}{e}_{\mathit{PYK}}\frac{{c}_{\mathit{adp}}}{{c}_{\mathit{atp}}}$ |
v _{ PDH } | ${k}_{\mathit{cat}}{e}_{\mathit{PDH}}\left(\frac{{c}_{\mathit{pyr}}^{4}}{{K}_{\mathit{pyr}}+{C}_{\mathit{pyr}}^{4}}\right)$ | v _{ PPC } | $\frac{{k}_{\mathit{cat}}{e}_{\mathit{PPC}}{c}_{\mathit{pep}}\left(1+\frac{{c}_{\mathit{fdp}}}{{K}_{\mathit{fdp}}}\right)}{{K}_{\mathit{pep}}+{C}_{\mathit{pep}}}$ |
v _{G 6PDH} | k _{ cat } e _{G 6PDH} c _{ nadp } c _{g 6p} | v _{ GND } | $\frac{{k}_{\mathit{cat}}{e}_{\mathit{GND}}{c}_{\mathit{nadp}}{c}_{6\mathit{pgc}}}{{c}_{\mathit{nadph}}{c}_{\mathit{atp}}\left({K}_{\mathit{pep}}+{C}_{\mathit{pep}}\right)}$ |
v _{ RPE } | ${k}_{\mathit{cat}}{e}_{\mathit{RPE}}\left({c}_{\mathit{ru}5\mathit{pD}}-\frac{{c}_{\mathit{xu}5\mathit{pD}}}{{K}_{\mathit{eq}}}\right)$ | v _{ RPI } | ${k}_{\mathit{cat}}{e}_{\mathit{RPI}}\left({c}_{\mathit{ru}5\mathit{pD}}-\frac{{c}_{r5p}}{{K}_{\mathit{eq}}}\right)$ |
v _{TKT 1} | ${k}_{\mathit{cat}}{e}_{\mathit{TKT}1}\left({c}_{r5p}{c}_{\mathit{xu}5\mathit{pD}}-\frac{{c}_{s7p}{c}_{\mathit{gap}}}{{K}_{\mathit{eq}}}\right)$ | v _{ TALA } | ${k}_{\mathit{cat}}{e}_{\mathit{TALA}}\left({c}_{\mathit{gap}}{c}_{s7p}-\frac{{c}_{e4p}{c}_{f6p}}{{K}_{\mathit{eq}}}\right)$ |
v _{TKT 2} | ${k}_{\mathit{cat}}{e}_{\mathit{TKT}2}\left({c}_{\mathit{xu}5\mathit{pD}}{c}_{e4p}-\frac{{c}_{f6p}{c}_{\mathit{gap}}}{{K}_{\mathit{eq}}}\right)$ |
Estimated parameter values and confidence intervals
Parameter | Value | 95% Confidence | Units | |
---|---|---|---|---|
PTS | k_{cat} | 20.7 | ±9.1 | mmol/mg protein/hr |
PGI | k_{cat} | 40.2 | ±9.5 | mmol/mg protein/hr |
K_{eq} | 1.23 | ±0.29 | Dimensionless | |
PFK | k_{cat} | 26 | ±35 | mmol/mg protein/mM/hr |
ALDO | k_{cat} | 3.965 | ±0.010 | mmol/mg protein/mM/hr |
K_{eq} | 0.18 | ±0.17 | mM | |
K_{fdp} | 0.0074 | ±0.0036 | mM | |
TPI | k_{cat} | 10000 | ±14000 | mmol/mg protein/hr |
K_{eq} | 0.11400 | ±0.00031 | Dimensionless | |
GAPD | k_{cat} | 10000 | ±4100 | mmol/mg protein/mM^{2}/hr |
K_{eq} | 1.21 | ±0.14 | Dimensionless | |
PGM | k_{cat} | 9995 | ±40 | mmol/mg protein/hr |
K_{eq} | 0.53570 | ±0.00063 | Dimensionless | |
PGK | k_{cat} | 54.3 | ±2.9 | mmol/mg protein/hr |
K_{eq} | 5512.1 | ±1.2 | Dimensionless | |
ENO | k_{cat} | 2 | ±45 | mmol/mg protein/mM/hr |
K_{eq} | 1.4 | ±4.1 | Dimensionless | |
PYK | k_{cat} | 40 | ±49 | mmol/mg protein/hr |
PDH | K_{cat} | 10.4 | ±4.6 | mmol/mg protein/hr |
K_{pyr} | 0.000020 | ±0.000029 | mM^{4} | |
PPC | k_{cat} | 2.15 | ±1.80 | mmol/mg protein/hr |
K_{fdp} | 2.5 | ±6.5 | mM | |
K_{pep} | 0.10 | ±0.16 | mM | |
G6PDH | k_{cat} | 859.6 | ±1.1 | mmol/mg protein/mM^{2}/hr |
GND | k_{cat} | 18.5 | ±10.8 | mmol*mM/mg protein/hr |
K_{6pg} | 0.021 | ±0.012 | mM | |
RPI | k_{cat} | 549.46 | ±0.69 | mmol/mg protein/mM/hr |
K_{eq} | 1.40000 | ±0.00019 | Dimensionless | |
RPE | k_{cat} | 10000 | ±13000 | mmol/mg protein/mM/hr |
K_{eq} | 0.4900 | ±0.0044 | Dimensionless | |
TKT1 | k_{cat} | 10000 | ±7800 | mmol/mg protein/mM^{2}/hr |
K_{eq} | 1.99 | ±0.012 | Dimensionless | |
TKT2 | k_{cat} | 10000 | ±5800 | mmol/mg protein/mM^{2}/hr |
K_{eq} | 3.500 | ±0.013 | Dimensionless | |
TALA | k_{cat} | 10000 | ±2300 | mmol/mg protein/mM^{2}/hr |
K_{eq} | 0.3675 | ±0.0021 | Dimensionless |
Three types of parameters remained in the simplified kinetic model: enzyme catalytic rate constants (k_{ cat }), equilibrium constants (K_{ eq }), and binding coefficients (K_{ m }). Equilibrium constants had the smallest relative confidence intervals, as changes in these parameters caused very large changes in the kinetic model enzyme and metabolite concentrations. Binding coefficients generally had small relative confidence intervals, as binding coefficients that could not be estimated reliably were eliminated in kinetic model simplification. Enzyme activities had the largest relative confidence intervals, and this was especially true for some reversible reactions. For reversible reactions that are often close to equilibrium (see next section), the catalytic rate constant is much larger than the reaction flux since the mechanistic function, f, is close to zero. The enzyme catalytic rate constants cannot be determined for these reactions with precision unless the reaction operates further from equilibrium. Disregarding the reactions that are often near equilibrium, the remaining enzyme activities have relative confidence intervals that are similar to binding and equilibrium constants.
Evaluation of optimal concentrations
Though parameters for these rate laws have been reported previously in the literature, we found that kinetic parameter values had to differ from these reported in vitro and in vivo values. When binding coefficient values were fixed to previously-reported in vitro or in vivo values (reported by Chassagnole et al.[17]) in the full kinetic model, no feasible solution could be found to Equations 1b-1e after several thousand nonlinear programming runs from random starting points where k_{ cat } values could vary. Biologically-relevant bounds on the equilibrium constants (K_{ eq }) and metabolite or protein concentrations (0.001 to 10 mM and 0.001 to 10 mg protein/gDCW, respectively) in part caused this infeasibility. When bounds were not included in the optimization problem, feasible solutions with large WSLS could be found. In contrast, feasible solutions with small WSLS were found after a few dozen starting points when binding coefficients and k_{ cat } values were all allowed to change simultaneously. This was not surprising, as there have been previous results showing that in vitro kinetic parameters are significantly different than their corresponding in vivo values [25-27].
The kinetic model predicted ribose-5-phosphate (r5p) and ribulose-5-phosphate (ru5p) concentrations deviated more from experimental measurements (Figure 3B) as compared to glycolytic and energy/redox carrier metabolites. This was not surprising because the pentose phosphate pathway metabolites had larger experimental variance, and thus were weighted less in WSLS parameter estimation. When all other measurements and kinetic rate laws were excluded from the WSLS problem, the rate laws for RPE and RPI were still difficult to fit to the r5p and ru5p measurements (data not shown), indicating that the proposed rate laws for these two reactions or their associated protein/metabolite measurements may be problematic.
Kinetic and thermodynamic limitations in central metabolism
Across all 20 considered conditions, the mean ΔG_{ ALDO }^{*} and ΔG_{ ENO }^{*} were less than −10 kJ/mol, indicating that fructose-bisphosphate aldolase (ALDO) and enolase (ENO) are often far from equilibrium and among the best targets to control flux through central metabolism [7]. The values for ΔG_{ ALDO }^{*} could be calculated directly from measurements as well, and these were also found to be on the same order of magnitude and far from equilibrium. These results indicate that, although ALDO and ENO are generally considered to be near equilibrium in human erythrocytes [29], this may change depending on the organism and environmental conditions. Other reactions, including pyruvate dehydrogenase (PDH), pyruvate kinase (PYK), phosphofructokinase (PFK), and glucose transport (PTS), are known to be irreversible and thus do not have a K_{ eq } in the kinetic model. These irreversible reactions were observed to have large negative ΔG_{ j }^{*} values in all experimental conditions when in vitro measured ΔG_{ j } values in standard conditions were used (data not shown) [21].
Most conditions had similar ΔG_{ j }^{*} values for a given reaction, indicating that the same thermodynamic limitations may arise regardless of the environmental or genetic changes made to the organism. However, five reactions (ribose phosphate isomerase, RPI; transketolase, TKT2; fructose bisphosphate aldolase, ALDO; phosphoglycerate kinase, PGK; and phosphoglucose isomerase, PGI) exhibited different ΔG_{ j }^{*} values for a few conditions, indicating that these five reactions may also control flux through central metabolism under certain conditions. However, the variation in estimated ΔG_{ j }^{*} values for the pentose phosphate reactions (RPI and TKT2) could be due to nosier (and/or missing) metabolite concentration measurements and to larger differences between predicted and measured concentrations for metabolites in this pathway (Figure 3). ΔG_{ j }^{*} variations could also be due to errors in estimated K_{ eq } values, since PGK and PGI both had large K_{ eq } confidence intervals (the confidence intervals for the other three reaction’s K_{ eq } values were small). For the conditions where the ΔG_{ j }^{*} values differed the most we did not observe any consistent patterns in the data that would explain these shifts in ΔG_{ j }^{*} values, such as high or low metabolite concentrations or high or low flux per enzyme concentration values.
Although many binding coefficient parameters were removed during kinetic model simplification, some binding coefficients could be estimated with reasonable confidence intervals. These binding coefficients could then be compared to metabolite concentrations, to identify enzymes whose activity appears to be substantially influenced by metabolite concentrations (Figure 5B). These results can be useful for metabolic engineering applications, as they show the enzyme-metabolite interactions that appear to be biologically-relevant in a variety of growth conditions. For example, inhibition of the glucose transport (PTS) by the pyruvate/phosphoenolpyruvate (c_{ pyr }/c_{ pep }) ratio was proposed in the rate laws by Chassagnole et al.[17] but this inhibition was unimportant and was removed during kinetic model simplification, while the proposed activation of phosphoenolpyruvate carboxylase (PPC) by fructose 1,6-bisphosphate (C_{ fdp }) was shown to substantially affect the PPC flux. Because of these kinetic model simplifications, the results clarify which inhibitory and activation effects are biologically-relevant in growing cells.
Predictions using constraint-based model with kinetic constraints
For the Δzwf, Δpgm, ΔtktB, and D=0.5 h^{-1} conditions, the mean residual from KFBA was smaller than the mean residual from FBA (Figure 6A). Uptake and secretion fluxes were also better predicted for these conditions. The Δpgm mutant has a growth defect (biomass yields are reduced by ~12% [16]), needing larger glucose uptake and carbon dioxide secretion rates when compared to the parental strain at the same growth rate. However, FBA predicts the same biomass yields for the Δpgm mutant as the parental strain at a growth rate of 0.2 h^{-1} and has 16.5% error and 27.7% error in the predicted glucose uptake and carbon dioxide production rates, respectively. The KFBA algorithm correctly predicts a lower Δpgm mutant biomass yield than the parental strain (Figure 6B), with just a 1% error and 3.2% error in the glucose uptake and carbon dioxide production rates, respectively. In this KFBA solution, PDH, PPC, and TKT2 reactions are all at their lower kinetic bounds, and carbon is forced through glycolysis and the pentose phosphate pathway that is not ultimately incorporated into biomass, and is instead secreted as carbon dioxide.
Not all predictions were improved by using the kinetic constraints. The residuals for the intracellular fluxes predicted by KFBA for the ΔgpmA mutant were greater than those predicted by FBA. In this case, the kinetic bounds associated with the PTS and PDH reactions caused KFBA to incorrectly predict fluxes through glycolysis, contributing to the large mean residual for KFBA. Removing the PTS and PDH bounds reduced the mean residual to 0.39, similar to FBA.
We also repeated the FBA and KFBA analysis using the iJO1366 genome-scale constraint-based model. The iJO1366 model contains more reactions in and around central metabolism. As a result, we found that fewer kinetic constraints need to be violated, as compared to the central model, and that pathways were predicted to be used that are not normally operational under glucose aerobic conditions (e.g., non-PTS glucose transport, glucose dehydrogenase, gluconate kinase, isocitrate lyase and xylose isomerase). We further evaluated how the KFBA solutions changed as the number of allowable kinetic flux bound violations (n*) increased (Figure 6C). These results show that there is a tradeoff between maximizing agreement with the kinetic flux bounds and not activating additional pathways (not included in the central metabolic model), which causes poorer agreement with experimental flux measurements. As observed with the central metabolic model, the biomass yields were predicted better when kinetic constraints were imposed (Figure 6D). To improve KFBA predictions of central metabolic fluxes, kinetic constraints for additional pathways (noted above) in the iJO1366 model need to be included.
Conclusions
Constraint-based modeling uses mass balances, flux capacities, reaction directionalities, and other relevant physical constraints to predict fluxes through metabolism. Although transcriptional regulatory and thermodynamic constraints have been integrated into this modeling approach, detailed kinetic constraints have not been extensively formulated, parameterized, and used in constraint-based models. Since kinetic constraints are often omitted from constraint-based models, some predicted flux distributions may not be achievable using native enzymes or protein levels. Incorporation of kinetic constraints into constraint-based allows multi-omic datasets to be used to find kinetic limitations on metabolic fluxes and suggests enzymes to target for improving cell behavior. For example, the Δpgm mutant is predicted to have lower biomass yields due to kinetic limitations, and the KFBA model suggests that decreasing PDH, PPC, and TKT2 levels would improve biomass yields for this mutant. One challenge with developing such kinetically-constrained models is finding kinetic parameter values that are consistent with experimental measurements.
In this study, a WSLS parameter estimation problem using multi-omic experimental data from a study by Ishii et al.[16] was formulated and solved. The parameter estimation results suggested changes to functional forms of rate laws, which were implemented to produce a simplified kinetic model. The simpler kinetic model is an improvement over the more detailed model because the parameters are better-resolved and the model could be solved more efficiently with a better fit to experimental data. Each of the retained metabolite-enzyme binding coefficient parameters were associated with measured metabolite concentrations, and more binding parameters could be retained in the future by measuring concentrations for more chemical species. Overall, the kinetic parameters we estimated could fit the kinetic model to 92.7% of the 720 measured metabolite and protein concentration measurements within one standard deviation across 20 different experiments.
The thermodynamic predictions about distance from equilibrium can also be used for metabolic engineering applications [7]. In this case, reactions that are far from equilibrium may be limited by kinetic regulation or enzyme levels, preventing them from reaching equilibrium. Here, we identified reactions that were far from equilibrium during growth in most conditions, which represent viable candidates for modifications to control flux through metabolism. The results of thermodynamic analysis from this study are consistent with other thermodynamic analysis of microbes. In a study by Kümmel et al.[7], it was found that the ALDO reaction in E. coli is not near equilibrium under physiological conditions. Klimacek et al.[30] found that the ALDO and ENO reactions are far from equilibrium in wild type and engineered Saccharomyces cerevisiae strains. Our results combined with these earlier studies indicate that E. coli and S. cerevisiae must use significantly different strategies for regulation of glycolysis than human erythrocytes, which instead closely regulate entry into and exit from glycolysis, and have reactions near equilibrium for all intermediate steps [29]. Enzymes were also identified where metabolite concentrations had significant binding, saturation, and allosteric regulatory effects. The conclusions about thermodynamics and enzyme binding are based on a simultaneous analysis of metabolomic, proteomic, and fluxomic data.
When kinetic constraints were imposed on a central metabolic constraint-based model the flux and biomass yield predictions were more accurately predicted by KFBA than FBA. Imposition of kinetic constraints in a genome-scale model provided mixed results, with more accurate biomass yields but worse overall flux residuals after incorporation of kinetic constraints. Since additional pathways are utilized in iJO1366 to match more kinetic flux bounds, the application of more kinetic constraints could improve predictions in more comprehensive models. The confidence intervals for kinetic parameters allowed reasonable flux ranges to be estimated from metabolite and enzyme concentration data. For the five test conditions that were evaluated, a median of 16 (out of 22) kinetic flux bounds were feasible. We considered potential errors in estimated parameter values to calculate the kinetic flux bounds, but errors in metabolite and protein concentration measurements could also be considered. The KFBA method used in this work could be applied to other reported rate laws, kinetic parameters, and concentration data, provided reasonable flux ranges can be calculated. However, identifying reasonable flux ranges may require parameter confidence intervals, which are not always available.
The presented approach to parameterize kinetic rate laws using in vivo data is general and can be applied to multi-omics data from other microbes. We chose published rate laws for a starting point for our kinetic model, but we note that these rate laws are almost entirely based on mass-action kinetics and Michaelis-Menten type inhibition. These rate laws were sufficient for our system, and we suggest the use of similar expressions for pathways where rate laws have not been proposed. The methods for imposing kinetic constraints in constraint-based models are also general, and can be used with rate laws with in vitro or in vivo determined parameters. Methods to parameterize and use kinetic rate laws in constraint-based models will benefit from more global and precise metabolomics and proteomics methods. Future work will involve including rate laws for other metabolic pathways (such as the citric acid cycle, fermentative and respiratory pathways) and estimating their in vivo kinetic parameters. Overall, this work illustrates how kinetic constraints can be used to improve predictions for intracellular fluxes and biomass yield and identify potential metabolic limitations through the integrated analysis of multi-omics datasets.
Author’s contributions
CC implemented the models and approach, performed the analysis, analyzed the data, and drafted the manuscript. JLR conceived of the study, participated in its design and coordination, and helped to analyze the data and draft the manuscript. All authors read and approved the final manuscript.
Declarations
Acknowledgements
This work was funded by the United States Department of Energy Great Lakes Bioenergy Research Center (DOE BER Office of Science DE-FC02-07ER64494). CC is also supported by a fellowship from the 3M Foundation. The authors also wish to acknowledge James B. Rawlings for helpful discussions.
Authors’ Affiliations
References
- Price ND, Reed JL, Palsson BO: Genome-scale models of microbial cells: evaluating the consequences of constraints. Nat Rev Microbiol 2004,2(11):886-897. 10.1038/nrmicro1023View ArticlePubMedGoogle Scholar
- Orth JD, Thiele I, Palsson BO: What is flux balance analysis? Nat Biotechnol 2010,28(3):245-248. 10.1038/nbt.1614PubMed CentralView ArticlePubMedGoogle Scholar
- Schuetz R, Kuepfer L, Sauer U: Systematic evaluation of objective functions for predicting intracellular fluxes in Escherichia coli. Mol Syst Biol 2007, 3: 119.PubMed CentralView ArticlePubMedGoogle Scholar
- Covert MW, Knight EM, Reed JL, Herrgard MJ, Palsson BO: Integrating high-throughput and computational data elucidates bacterial networks. Nature 2004,429(6987):92-96. 10.1038/nature02456View ArticlePubMedGoogle Scholar
- Chandrasekaran S, Price ND: Probabilistic integrative modeling of genome-scale metabolic and regulatory networks in Escherichia coli and Mycobacterium tuberculosis. Proc Natl Acad Sci USA 2010,107(41):17845-17850. 10.1073/pnas.1005139107PubMed CentralView ArticlePubMedGoogle Scholar
- Shlomi T, Eisenberg Y, Sharan R, Ruppin E: A genome-scale computational study of the interplay between transcriptional regulation and metabolism. Mol Syst Biol 2007, 3: 101.PubMed CentralView ArticlePubMedGoogle Scholar
- Kümmel A, Panke S, Heinemann M: Putative regulatory sites unraveled by network-embedded thermodynamic analysis of metabolome data. Mol Syst Biol 2006, 2: 2006-0034.PubMed CentralView ArticlePubMedGoogle Scholar
- Henry CS, Broadbelt LJ, Hatzimanikatis V: Thermodynamics-based metabolic flux analysis. Biophys J 2007,92(5):1792-1805. 10.1529/biophysj.106.093138PubMed CentralView ArticlePubMedGoogle Scholar
- Hoppe A, Hoffmann S, Holzhutter HG: Including metabolite concentrations into flux balance analysis: thermodynamic realizability as a constraint on flux distributions in metabolic networks. BMC Syst Biol 2007, 1: 23. 10.1186/1752-0509-1-23PubMed CentralView ArticlePubMedGoogle Scholar
- Henry CS, Jankowski MD, Broadbelt LJ, Hatzimanikatis V: Genome-scale thermodynamic analysis of Escherichia coli metabolism. Biophys J 2006,90(4):1453-1461. 10.1529/biophysj.105.071720PubMed CentralView ArticlePubMedGoogle Scholar
- Beg QK, Vazquez A, Ernst J, de Menezes MA, Bar-Joseph Z, Barabasi AL, Oltvai ZN: Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity. Proc Natl Acad Sci USA 2007,104(31):12663-12668. 10.1073/pnas.0609845104PubMed CentralView ArticlePubMedGoogle Scholar
- Yizhak K, Benyamini T, Liebermeister W, Ruppin E, Shlomi T: Integrating quantitative proteomics and metabolomics with a genome-scale metabolic network model. Bioinformatics 2010,26(12):i255-260. 10.1093/bioinformatics/btq183PubMed CentralView ArticlePubMedGoogle Scholar
- Hanly TJ, Henson MA: Dynamic flux balance modeling of microbial co-cultures for efficient batch fermentation of glucose and xylose mixtures. Biotechnol Bioeng 2010,108(2):376-385.View ArticleGoogle Scholar
- Mahadevan R, Edwards JS, Doyle FJ 3rd: Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J 2002,83(3):1331-1340. 10.1016/S0006-3495(02)73903-9PubMed CentralView ArticlePubMedGoogle Scholar
- Zhuang K, Izallalen M, Mouser P, Richter H, Risso C, Mahadevan R, Lovley DR: Genome-scale dynamic modeling of the competition between Rhodoferax and Geobacter in anoxic subsurface environments. ISME J 2010,5(2):305-316.PubMed CentralView ArticlePubMedGoogle Scholar
- Ishii N, Nakahigashi K, Baba T, Robert M, Soga T, Kanai A, Hirasawa T, Naba M, Hirai K, Hoque A, et al.: Multiple high-throughput analyses monitor the response of E. coli to perturbations. Science 2007,316(5824):593-597. 10.1126/science.1132067View ArticlePubMedGoogle Scholar
- Chassagnole C, Noisommit-Rizzi N, Schmid JW, Mauch K, Reuss M: Dynamic modeling of the central carbon metabolism of Escherichia coli. Biotechnol Bioeng 2002,79(1):53-73. 10.1002/bit.10288View ArticlePubMedGoogle Scholar
- Palsson B: Systems biology: properties of reconstructed networks. Cambridge; New York: Cambridge University Press; 2006.View ArticleGoogle Scholar
- Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, Palsson BO: A comprehensive genome-scale reconstruction of Escherichia coli metabolism-2011. Mol Syst Biol 2011, 7: 535.PubMed CentralView ArticlePubMedGoogle Scholar
- Chapman AG, Fall L, Atkinson DE: Adenylate energy charge in Escherichia coli during growth and starvation. J Bacteriol 1971,108(3):1072-1086.PubMed CentralPubMedGoogle Scholar
- Goldberg RN, Tewari YB, Bhat TN: Thermodynamics of enzyme-catalyzed reactions-a database for quantitative biochemistry. Bioinformatics 2004,20(16):2874-2877. 10.1093/bioinformatics/bth314View ArticlePubMedGoogle Scholar
- Antoniewicz MR, Kelleher JK, Stephanopoulos G: Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements. Metab Eng 2006,8(4):324-337. 10.1016/j.ymben.2006.01.004View ArticlePubMedGoogle Scholar
- Lewis NE, Hixson KK, Conrad TM, Lerman JA, Charusanti P, Polpitiya AD, Adkins JN, Schramm G, Purvine SO, Lopez-Ferrer D, et al.: Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol 2010, 6: 390.PubMed CentralView ArticlePubMedGoogle Scholar
- Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally sloppy parameter sensitivities in systems biology models. PLoS Comput Biol 2007,3(10):1871-1878.View ArticlePubMedGoogle Scholar
- Van Eunen K, Kiewiet JA, Westerhoff HV, Bakker BM: Testing biochemistry revisited: how in vivo metabolism can be understood from in vitro enzyme kinetics. PLoS Comput Biol 2012,8(4):e1002483. 10.1371/journal.pcbi.1002483PubMed CentralView ArticlePubMedGoogle Scholar
- Huang X, Holden HM, Raushel FM: Channeling of substrates and intermediates in enzyme-catalyzed reactions. Annu Rev Biochem 2001, 70: 149-180. 10.1146/annurev.biochem.70.1.149View ArticlePubMedGoogle Scholar
- Ellis RJ: Macromolecular crowding: obvious but underappreciated. Trends Biochem Sci 2001,26(10):597-604. 10.1016/S0968-0004(01)01938-7View ArticlePubMedGoogle Scholar
- Wang L, Birol I, Hatzimanikatis V: Metabolic control analysis under uncertainty: framework development and case studies. Biophys J 2004,87(6):3750-3763. 10.1529/biophysj.104.048090PubMed CentralView ArticlePubMedGoogle Scholar
- Garrett R, Grisham CM: Biochemistry. Fort Worth: Saunders College Pub; 1995.Google Scholar
- Klimacek M, Krahulec S, Sauer U, Nidetzky B: Limitations in xylose-fermenting Saccharomyces cerevisiae, made evident through comprehensive metabolite profiling and thermodynamic analysis. Appl Environ Microbiol 2010,76(22):7566-7574. 10.1128/AEM.01787-10PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.