Bridging the gap between gene expression and metabolic phenotype via kinetic models
© Vital-Lopez et al.; licensee BioMed Central Ltd. 2013
Received: 27 September 2012
Accepted: 27 June 2013
Published: 22 July 2013
Despite the close association between gene expression and metabolism, experimental evidence shows that gene expression levels alone cannot predict metabolic phenotypes, indicating a knowledge gap in our understanding of how these processes are connected. Here, we present a method that integrates transcriptome, fluxome, and metabolome data using kinetic models to create a mechanistic link between gene expression and metabolism.
We developed a modeling framework to construct kinetic models that connect the transcriptional and metabolic responses of a cell to exogenous perturbations. The framework allowed us to avoid extensive experimental characterization, literature mining, and optimization problems by estimating most model parameters directly from fluxome and transcriptome data. We applied the framework to investigate how gene expression changes led to observed phenotypic alterations of Saccharomyces cerevisiae treated with weak organic acids (i.e., acetate, benzoate, propionate, or sorbate) and the histidine synthesis inhibitor 3-aminotriazole under steady-state conditions. We found that the transcriptional response led to alterations in yeast metabolism that mimicked measured metabolic fluxes and concentration changes. Further analyses generated mechanistic insights of how S. cerevisiae responds to these stresses. In particular, these results suggest that S. cerevisiae uses different regulation strategies for responding to these insults: regulation of two reactions accounted for most of the tolerance to the four weak organic acids, whereas the response to 3-aminotriazole was distributed among multiple reactions. Moreover, we observed that the magnitude of the gene expression changes was not directly correlated with their effect on the ability of S. cerevisiae to grow under these treatments. In addition, we identified another potential mechanism of action of 3-aminotriazole associated with the depletion of tetrahydrofolate.
Our simulation results show that the modeling framework provided an accurate mechanistic link between gene expression and cellular metabolism. The proposed method allowed us to integrate transcriptome, fluxome, and metabolome data to determine and interpret important features of the physiological response of yeast to stresses. Importantly, given its flexibility and robustness, our approach can be applied to investigate the transcriptional-metabolic response in other cellular systems of medical and industrial relevance.
KeywordsGene expression Kinetic models Metabolic networks S. cerevisiae Transcriptomics Fluxomics Metabolomics
It is well known that cells regulate gene expression to carry out different functions depending on their physiological state and environment. However, it is less well understood how this regulation is orchestrated and how gene expression changes drive cells to adapt particular phenotypes. Developments in high-throughput technologies have contributed to answer these questions by generating a wealth of data on different cellular components and processes (e.g., transcriptome, proteome, metabolome, fluxome, and protein-protein interaction data). Hence, one of the challenges in systems biology is how to integrate and analyze such data to elucidate the underlying cellular physiology. In particular, the development of genome-scale computational models and analysis tools can help expand our understanding of how gene transcription alters cellular metabolism.
Different approaches have already made considerable headway in integrating gene expression and metabolism [1–6]. Perhaps the most developed efforts are based on combining stoichiometric models of metabolic networks and gene expression data. In these approaches, gene expression levels are used to parameterize the flux capacity of metabolic reactions to create context-specific models [7–9]. For example, we followed this approach to characterize the metabolic adaptations of Mycobacterium tuberculosis to hypoxia and identify metabolic alterations required for adaptation to a lifestyle of low metabolic activity . Alternatively, computational approaches have been developed to infer regulatory networks from gene expression data , which in turn have been integrated with metabolic network models to describe the adaptation of an organism to different conditions [12–15].
Combining stoichiometric models of metabolic networks and gene expression data has proven useful in analyzing transcriptome, proteome, and fluxome data but presents limitations in analyzing metabolome data. These limitations can be overcome using kinetic models, in which metabolite concentrations are the primary variables as opposed to fluxes in constraint-based methods. However, the use of large-scale kinetic models (i.e., with hundreds of reactions) has been daunted by the general belief that the chances of obtaining a useful model, given the lack of accurate reaction rate expressions and kinetic parameters, are low. This paradigm has begun to change due, in part, to the high-throughput techniques that have increased the abundance, quality, and scope of the data needed for model construction.
In addition to data availability, there are two other factors, arising from the biology of the systems, that ease the construction of large-scale kinetic models . The first one is the observation that the structure of a biological network (i.e., what the network components are and how they are connected) largely determines its function, as observed in constraint-based analyses . Thus, the available reconstructions of metabolic networks provide us with more than a solid scaffold to construct kinetic models: the performance of the network is confined within well-characterized limits. The second factor is the “sloppiness” of parameter sensitivities, which seems to be a widespread property of models of biological systems . This sloppiness property implies that most of the model parameters cannot be collectively estimated with certainty, even by fitting large amounts of “ideal” data. Paradoxically, it also implies that knowledge of the precise value of most parameters is not critical for describing a system’s behavior. Motivated by these factors, methods to construct large-scale kinetic models of metabolism have started to emerge [19–22].
In this work, our objective was to investigate how the response of a cell to a perturbation (in terms of transcriptome or proteome data) induces changes in its phenotype (in terms of fluxome and metabolome data). For this purpose, we developed a computational approach based on kinetic models that provides a mechanistic link between transcriptional regulation and metabolism. Our proposed modeling framework overcomes the major obstacles in the construction of large-scale kinetic models of metabolism, namely, the detailed definition of appropriate reaction rate expressions and the determination of model parameters. As in previous approaches [19–22], we automatically translated a metabolic network model into a kinetic model using generic expressions, a particular case of generalized mass action (GMA) kinetics, for the reaction rates . However, in contrast to these approaches, our method does not require extensive parameter estimation, mining the literature, or using random-sampling schemes to obtain parameter values. Most of the model parameters are obtained directly from experimental data that are routinely available (i.e., protein or gene expression data and flux distributions or uptake/production rates). Although the models could be used to investigate dynamic behavior, this would require additional input parameters in terms of an extensive set of metabolite concentrations. However, as these data are typically not available, and similar to other approaches such as ensemble modeling , we have used the proposed models to describe and analyze steady-state behavior.
Here, we constructed kinetic models to analyze the steady-state metabolism of S. cerevisiae based on two independent studies in which the transcriptional and metabolic responses to treatment were measured in chemostat cultures with weak organic acids (WOAs)  and under histidine starvation conditions . The simulation results demonstrated that integration of gene expression with metabolic network models enabled us to capture aspects of the response of S. cerevisiae that would not have been possible by the independent analyses of the gene expression data or the metabolic network alone.
where β is a parameter that relates the rate of the forward and backward reactions to the overall flux at the reference condition. The value of β depends on the equilibrium constant and on the reactants concentration at the reference condition. If these data are not available, as in the experiments analyzed here, β could be estimated by a fitting procedure using other available data. However, to avoid over fitting and assuming that the model behavior has low sensitivity to individual β s we started by approximating all β s with the same value (except for reactions in parallel routes that must satisfy additional thermodynamic constraints as detailed in Additional file 1) and checked whether it was necessary to estimate individual β s. As shown in the Results Section, there was no need for estimating individual β s as a single-value approximation proved satisfactory.
where μ denotes the growth rate, λ denotes a correction factor, X MW denotes the biomass molecular weight, ϕ i denotes the moles of carbon per mole of biomass precursor, and the summation included only the drain fluxes to biomass. Note that this definition of the biomass growth rate can predict a non-zero rate even if some of the drain fluxes are zero. However, such extreme cases were not observed in the simulations carried out in this paper.
where the subscript i comprises the metabolites consumed in each drain flux.
where C is a diagonal matrix with elements equal to the absolute metabolite concentrations used for normalization, c represents the vector of normalized metabolite concentrations and ċ denotes its time derivative, S denotes the stoichiometric matrix of the metabolic network reconstruction, r represents the vector of reaction rates, v denotes the flux distribution, g represents the vector of gene expression ratios, and p denotes a vector of the other model parameters (i.e., α, β, λ, and other condition-specific parameters). Under steady-state conditions, C is not required and, thus, for steady-state analysis, the only parameters to be estimated are v, g, and p.
In Step 2, we parameterized the model for the reference condition. Using the reference condition for normalizing the metabolite concentrations and gene expression levels, both c and g become equal to 1.0, and r=v ref , where v ref is the flux distribution at the reference condition. Therefore, for steady-state analysis, the model for the reference condition was parameterized with v ref . A flux distribution determined using 13C-labeling experiments provides a good estimate of v ref . If such flux distribution is not available, a reasonable estimate can be obtained using exchange fluxes, as described in Additional file 1.
Another notable feature of the method is that the model can be parameterized to simulate other conditions using the gene expression ratio between the condition of interest and the reference condition (Step 3). We assumed that relative changes in gene expression led to similar relative changes in protein abundance and we neglected post-translational and other regulatory mechanisms of enzymatic activity. Note that, if available, proteome data can be used instead of gene expression data. For reactions associated with multiple genes, we computed an overall gene expression change as described in Additional file 1.
In Step 4, we tuned the constructed models by comparing model predictions (i.e., metabolic fluxes and metabolite concentration changes) with experimental measurements. We resolved detected inconsistencies until the model gave satisfactory performance. Finally, once its performance is satisfactory, the model can be used to carry out different model-based analyses, such as predicting non-measured variables, determining the effect of a particular expression level on a given metabolic function, or to identify important reactions in the network (Step 5).
The experimental data for the analysis of the response of S. cerevisiae to treatment with a WOA (acetic acid, benzoic acid, propionic acid, or sorbic acid) were obtained from Abbott et al. . To compute the gene expression ratios from the raw intensity values, the microarray data were scaled such that the average intensity for each microarray was 150.0. For each condition, the median intensity (of three replicates) was used as the expression level of every gene in each condition. Every treatment-reference condition pair was smoothed using the Lowess smoothing method  and the gene expression ratios were computed from the smoothed data. The uptake rate of glucose and production rates of ethanol, glycerol, lactic acid, acetic acid, and CO2 were reported, as well as the biomass yield and concentration, and the extracellular glucose concentration for each treatment condition.
For the analysis of the metabolic response induced by the Gcn4 regulator, the experimental data were taken from Moxley et al. . They report the gene expression data, flux distribution estimated from 13C-labeling experiments, and the concentration changes of 17 free amino acids for chemostat culture of wild-type and gcn4-knockout mutant strains of S. cerevisiae.
Implementation and availability
Model construction, data processing, and simulations were carried out in MATLAB (2011b, The MathWorks Inc., Natick, MA). The kinetic model (MATLAB scripts and in SBML format) and parameter sets for simulating both experiments are provided in Additional file 2.
Construction of large-scale kinetic models
Fitting parameters for the kinetic model of S. cerevisiae metabolic network
Kinetic order for biomass drain fluxes
Ratio of the forward (or backward if reference flux is negative) reaction rate
to the overall rate of the reversible reactions
Correction factor for biomass growth rate
Inhibition factor of the synthesis of glycine from serine by 3-AT
Inhibition factor of the synthesis of histidine 3-AT
WOA uptake rate
Response of S. cerevisiae to histidine starvation
The activator protein Gcn4 of S. cerevisiae regulates the expression of nearly all genes encoding enzymes involved in amino acid synthesis under starvation conditions . Moxley et al.  studied the regulatory and metabolic changes induced by Gcn4 under histidine-deficient conditions. Specifically, they cultivated wild-type and gcn4-knockout mutant (Δgcn4) strains of S. cerevisiae in aerobic chemostats treated with 3-aminotriazole (3-AT), an inhibitor of imidazoleglycerol-phosphate dehydratase, the sixth step of the histidine synthesis pathway. The concentration of 3-AT was adjusted such that the Δgcn4 and wild-type cultures produced similar biomass levels and uptake and production rates of extracellular metabolites. They measured gene expression levels (See Additional file 1 for details about gene expression data processing), metabolic fluxes (using 13C-labeling experiments), and the intracellular concentration of free amino acids for each culture. We used the measured metabolic fluxes of the Δgcn4 culture as the reference condition v ref and the gene expression ratios between the two cultures g to parameterize the model for simulating the wild-type cultures. The parameters in p used in these simulations are given in Table 1. First, we showed that the model was able to reproduce the experimental data. We then examined possible mechanisms of action of 3-AT and delineated the effect of the gene expression changes on yeast’s ability to grow when exposed to this drug.
The model quantitatively links gene expression regulation with metabolism
Note that although we used the experimental concentration changes to estimate the value of the fitting parameters, the high ρ of 0.96 could not be achieved without parameterizing the model with the gene expression data. Nonetheless, the ρ without using the gene expression data was relatively high (ρ = 0.80; see Additional file 1). Based on these results, we can derive two general observations. First, the structure of the metabolic network, which we exploited to constrain the kinetic parameters v in our modeling framework, considerably contributed to the explanation of the experimental observations as we initially assumed. Second, gene expression changes were required to further improve the simulation results. Thus, both the metabolic network and the gene expression changes were required by the framework to establish the mechanistic link between gene expression regulation and metabolism.
Model-based identification of mechanisms of action of 3-AT
Predicted targets of 3-aminotriazole in the S. cerevisiae central carbon metabolic network
SER ↔ GLY + METTHF (forward)
SER ↔ GLY + METTHF (backward)
R5P + METTHF + 2 ATP → IAP
HIS + ATP → HISBIO
IAP + GLU → HIS + AKG + 2 NADH
G3P ↔ PEP + NADH + ATP (forward)
G3P ↔ PEP + NADH + ATP (backward)
NADH → ATP
HSER + METTHF + ACCOA + 2 ATP → MET + ACE + 3 NADH
MET + ATP → METBIO
Effect of gene expression regulation on the tolerance to 3-AT treatment
Another interesting result was that the magnitude of the individual effects was not correlated with the magnitude of the gene expression changes (ρ = 0.06). Moreover, nine of the top ten reactions in Figure 4 had associated gene expression changes of less than two-fold. This suggests that the magnitude of gene expression changes may be a poor predictor of their importance, supporting the notion that analyses biased towards large gene expression changes may miss important insights. Note, however, that in general, small gene expression changes have more uncertainty and are more sensitive to normalization errors than large expression changes.
Modeling the response of S. cerevisiae to treatment with WOAs
Constructed models captured S. cerevisiae response to WOA treatment
Predicted metabolic responses to weak organic acid treatments
Effect of gene expression data (GED) on the predicted concentrations under weak organic acid treatments
We further investigated the effect of the transcriptional response by comparing the predicted metabolic response of untreated cultures while considering the gene expression levels of the untreated or the treated cultures (i.e., simulations without gene expression data vs. simulations with gene expression data). Under the reference condition, the WOA uptake rate is the corresponding diffusive uptake flux of the acetic acid produced by S. cerevisiae. Figure 6B shows the normalized SSE of the predicted exchange fluxes and biomass yield with respect to the experimental values for the treatment conditions. In contrast to simulations without gene expression data, which by construction simulated the reference condition and had a normalized SSE of 1.0, predictions using gene expression were closer to the experimental values of the treatment conditions. This shows that the gene expression data captured, to a certain degree, the metabolic response of S. cerevisiae to the WOA treatments. Taken together, these results suggest that while the network structure had a predominant role on the metabolic flux distribution, gene expression changes contributed to flux regulation and had a major effect on metabolite concentration changes.
Effect of gene expression changes on WOA tolerance
Identification of key gene expression changes for tolerance to WOA
Contribution to tolerance to weak organic acid (WOA) treatment and gene expression changes of the two most important reactions
EC 126.96.36.199 (Glucose → Glucose-6P)
EC 188.8.131.52 (Pyruvate → Acetaldehyde)
Sensitivity and robustness of model predictions
An advantage of our method is that it has a small number of general fitting parameters: α, β, and γ. We investigated the effect of these fitting parameters on model predictions by simulating the model for different values of these parameters. For the WOA treatment experiments, we determined the normalized SSE between the predicted and experimental values of the exchange fluxes and biomass yield as a function of each parameter. The analysis showed that the predicted fluxes were robust with respect to changes in these parameters (see Additional File 1). The simulations of the histidine starvation experiments were relatively robust to variations in these parameters around the values reported in Table 1, as shown in Additional File 1.
We also investigated if the method was sensitive to the input gene expression data or if the results could be obtained with arbitrary data. Briefly, we simulated the models with data generated by randomly shuffling the original expression data. The simulation results showed that it is unlikely that similar results could be obtained using random gene expression data (see Additional File 1). Moreover, uncertainty propagation analysis showed that the method is robust with respect to experimental noise in the gene expression data and flux distributions (Additional File 1).
Long-standing barriers impeding the construction of large-scale kinetic models of metabolism are being overcome with the help of developments in high-throughput technologies and computational analyses. Modelers are now faced with the challenge of integrating the increasingly available building blocks to create coherent mathematical representations of biological systems. Here, we presented our efforts to develop a modeling framework for constructing large-scale kinetic models that mechanistically link transcriptional regulation and metabolism. This allowed us to gain understanding of complex physiological relations from fluxome, metabolome, and gene expression data. We demonstrated the ability of our method to capture these relations, its flexibility to simulate different experiments, and its robustness with respect to modeling approximations and data uncertainty by analyzing the response of S. cerevisiae under different stress conditions. Importantly, our approach can be applied to other organisms of medical and industrial relevance (or cell types in multi-cellular organisms) for which a metabolic network reconstruction, metabolic flux measurements, and gene expression data are available for the conditions of interest.
The method provides efficient solutions to large-scale modeling challenges
One of the major challenges in constructing large-scale kinetic models is the definition of appropriate reaction rate expressions. Instead of defining mechanistic reaction rate expressions on a case-by-case basis, some approaches streamline this process by relying on generic expressions to translate a metabolic network into a kinetic model in an automated or semi-automated fashion. Different general forms have been proposed, such as log-linear kinetics , Michaelis-Menten-type kinetics , “convenience” kinetics , or GMA kinetics . GMA kinetics are used, for example, in ensemble modeling  and mass action stoichiometric simulation (MASS) models . In ensemble modeling and MASS models, the enzymatic reactions are decomposed into their elementary steps, and each step is then modeled using mass action kinetics. The decomposition increases the resolution of the model, preserves enzyme saturation behavior, and simplifies the parameter estimation problem, but at the price of considerably increasing the size of the model (i.e., the number of dynamic variables and model parameters) and the amount of data required to estimate parameter values. In contrast, we used a special case of GMA kinetics that requires a minimal number of parameters, which can be obtained directly from available experimental data (see Methods and Additional File 1). Moreover, enzymatic reactions were not decomposed into elementary steps to avoid increasing the size of the model.
Another challenge is the determination of model parameter values. The difficulty in solving this problem is linked to the form of the kinetic expressions and to the availability of experimental data. If experimental data are not available, approaches such as log-linear kinetics and “convenience” kinetics require mining the literature for parameter values, which (aside from the skepticism about the validity of combining parameter values from different conditions to simulate a specific experiment) could be impractical for large-scale models. Approaches using GMA kinetics partially avoid literature mining. In these approaches, such as MASS modeling , thermodynamic information collected from the literature (e.g., equilibrium constants, Gibbs free energies, etc.) is combined with experimentally determined metabolite and/or enzyme concentrations and flux distributions to estimate the remaining model parameters (i.e., the rate constants). For the common case of incomplete data, the missing information is approximated to “typical” values or is randomly generated to create an ensemble of models that are screened for models that agree with experimental observations . Based on the “sloppiness” property, we would expect that models parameterized using “typical” values will perform reasonable well. However, the typical values generally fall within relatively wide ranges, making the selection of parameter values to simulate a particular condition a non-trivial task. In contrast, the rate expressions we used enabled us to readily obtain the bulk of the model parameters (221 out of 227) directly from available experimental data (i.e., flux distributions and gene expression; see Methods and Additional File 1). Moreover, we circumvented mining the literature or using randomly generated values for thermodynamic parameters by assuming a single parameter (β) for relating the forward and backward reaction rates to the overall rate for all reversible reactions. This crude approximation, inspired in part by the “sloppiness” property of biological systems, worked surprisingly well for the examples studied here. Our method performed well even if the uptake and production rates of extracellular metabolites were the only metabolic data available, as demonstrated in the analysis of S. cerevisiae tolerance to WOAs (see Additional File 1 for simulations using only uptake and production rates of extracellular metabolites under histidine starvation).
An additional attribute of our method is the use of gene expression data to parameterize the model to simulate different conditions, an element that has been used in constraint-based approaches to create context-specific models [7–10], but which has not been fully exploited in other kinetic modeling approaches. An exception is the work by Bruck et al. , in which gene expression was integrated with a kinetic model of S. cerevisiae glycolysis based on a mechanistic model developed by Teusink et al. . However, Bruck et al.  estimated a subset of 31 parameters to fit the model to data from all conditions they simulated and did not present simulations without the gene expression data, preventing an assessment of the contribution of gene expression changes. In contrast, our models were able to simulate metabolic responses with a smaller subset of fitting parameters and our analysis showed the important role of gene expression on model predictions. Note that requiring gene expression data in order to simulate other conditions could also be considered a weakness, but no other model includes the prediction of protein/gene expression changes for the systems of the size of the network we analyzed.
Constructed models generated biological insights
We demonstrated that the constructed models were able to integrate transcriptional and metabolic responses to produce insights that would have been difficult to grasp from the analysis of the individual responses. For example, in their analysis of S. cerevisiae response to WOAs, Abbott et al.  identified differentially expressed genes as those with an expression change larger than two-fold and a false discovery rate lower than 0.5%. With these criteria, they found hundreds of differentially expressed genes under each treatment condition, but only 14 genes that were upregulated under all treatment conditions. Therefore, they concluded that the generic (i.e., common to all treatments) transcriptional response to WOAs was minimal and suggested that more relevance should be given to the specific responses to the specific treatment conditions. We agree that attention should be paid to the specific responses, but our analysis also suggests that the generic response, despite involving a few genes, is a major factor contributing to WOA tolerance. Based on our simulation results, we hypothesize that S. cerevisiae tightly regulates the expression levels of two reactions (glucose uptake-phosphorylation and decarboxylation of pyruvate to acetaldehyde) to increase the tolerance under all treatment conditions. Firstly, this generic response was not identified in the Abbott et al.  analysis because the gene expression changes for these reactions did not meet their criteria for differentially expressed genes (see Additional File 1). Secondly, we estimated that regulating these two reactions accounted for most of the increase in tolerance to WOAs (Figure 8 and Table 5). If correct, this hypothesis implies that S. cerevisiae has a generic response to WOAs that is critical for the adaptation to these stressors.
Identification of important reactions in a metabolic network has been one of the major goals of several model-based approaches. For example, Kummel et al.  developed a thermodynamics-based method to identify regulated reactions, assuming that reactions far from equilibrium are more likely to be regulated. In contrast with our method, their approach does not use any kinetic information but requires thermodynamic and metabolome data. In another example, Smallbone et al.  combined log-linear kinetics with metabolic control analysis [37, 38] to identify reactions exerting the most control over biomass production in a genome-scale metabolic network of S. cerevisiae. Similar to these efforts, our method was able to identify important regulated reactions under specific conditions. However, our method also provided mechanistic insights into how the cell regulates such reactions through transcriptional regulation and how this response is reflected in its phenotype.
In another effort to link the regulatory and metabolic responses, Moxley et al.  proposed a hybrid approach to predict changes in metabolic fluxes using gene expression changes. Their approach was based on the assumption that gene expression changes and fluxes are more correlated in pathways with fewer metabolite-enzyme interactions (metabolite-enzyme interactions exist between an enzyme and metabolites that regulate its activity). Thus, their approach combined a metabolic network model with a metabolite-enzyme interaction network. Using this approach, they predicted flux changes that had a relatively high correlation (ρ = 0.80) with the experimentally estimated flux changes for a subset of reactions. For the same subset, our model predictions showed a considerably higher correlation (ρ = 0.96). Moreover, our method required less information because knowledge of the metabolite-enzyme interaction network is not needed. Interestingly, their predictions, using only the metabolic network model (without considering metabolite-enzyme interactions), had a similar ρ of approximately 0.75, reflecting the major contribution of the network structure to its function. In terms of biological insights, they observed a redistribution of the glycine synthesis fluxes. They proposed that the increase in glycine production from threonine is mediated by the increased expression of the associated genes, but they do not fully explain why the flux from serine to glycine decreased. Our analysis led to the plausible explanation that the decrease in the flux from serine to glycine could have been caused by the decrease of tetrahydrofolate, which, in turn, could have been caused by off-target inhibitions of 3-AT. In addition, and in contrast with their approach, our method also predicted concentration changes. In fact, we are unaware of other modeling efforts with similar scope that produce similar levels of accuracy, using condition-specific data directly as model parameters and using only five fitting parameters.
An additional conjecture about the use of gene expression changes to parameterize protein activity changes can be derived from our simulation results. We omitted post-translational and other regulatory mechanisms and yet the model predictions were consistent with experimental data. This suggests that, for the metabolic network and the experiments considered here, transcriptional regulation was the main mechanism that regulated the response at the system level. Moreover, the accuracy of the model predictions suggests that gene expression changes were a good approximation for protein level changes, in agreement with experimental observations [27, 39].
The proposed method does not need knowledge of the absolute values of metabolite concentrations for steady-state simulations, but these are required for analysis of transient behavior. Developments in analytical techniques have increased the accuracy and scope of metabolite concentration measurements. However, such data are still generally incomplete and, thus, missing data must be estimated or assumed. Note that the requirement of metabolite concentrations to describe dynamic behavior is common to similar modeling approaches. Thus, it remains to be investigated how the proposed modeling framework performs in describing dynamic and transient properties associated with metabolic processes.
The models constructed with the proposed method present some limitations. For example, the generic rate expressions may be poor approximations for some reactions or may miss important allosteric regulations (e.g., feedback loops) and other factors that have an effect on protein activity and abundance (e.g., post-translational modifications). Lumping sequential reactions reduced the size of the model. However, in our approach, the rate expressions for lumped reactions are only an approximation to the sequence of individual reactions. In the experiments we analyzed, the final results were not sensitive to our somewhat arbitrary parameter choice m i and β. This may not be always the case and estimating more accurate parameters values may be necessary. As for any method, identifying and correcting modeling errors is a painstaking task. This could be especially true for automated model generation. Procedures to address this problem in a systematic way need to be developed. Furthermore, our method needs to be tested to determine whether it can be applied to genome-scale metabolic networks. Such application could be problematic because of the higher uncertainty of lowly expressed genes and small metabolic fluxes, the buildup of approximation errors, and numerical challenges to solve the model. Regarding its scope, the proposed method is limited to gene expression and metabolism. Although it enables a deeper, mechanistic analysis of these processes, further developments to include other cellular processes (e.g., signal transduction, cell division, etc.) would greatly enhance the modeling framework.
In summary, we investigated how gene expression changes induce metabolic responses when cells adapt to a stressful condition. For this purpose, we developed a modeling framework for constructing and simulating large-scale kinetic models that provided a mechanistic link between transcriptional regulation and cellular metabolism. Analysis of the response of S. cerevisiae to treatment with WOA and under histidine starvation generated several insights and testable hypotheses: 1) 3-AT also inhibits the synthesis of tetrahydrofolate; 2) S. cerevisiae has an important generic response to WOA, involving glucose uptake and decarboxylation of pyruvate to acetaldehyde; 3) the contribution to tolerance to 3-AT is distributed among several reactions while the contribution to tolerance to WOA is mainly concentrated in two reactions; 4) the magnitude of gene expression changes was not correlated with the magnitude of their effect on the overall response. Taken together, these results suggest that the proposed framework is able to dissect different “omics” data to determine important features of the transcriptional-metabolic response of S. cerevisiae.
All authors contributed to the design of the research. FGV performed the computational implementations. FGV and AW wrote the article, which was edited by JR. All authors read and approved the final manuscript.
The authors thank Dr. M. D. AbdulHameed for his help in identifying possible targets of 3-AT, and Dr. X. Fang and Dr. B. Dutta for valuable discussions.
- Feist AM, Herrgard MJ, Thiele I, Reed JL, Palsson BO: Reconstruction of biochemical networks in microorganisms. Nat Rev Microbiol. 2009, 7: 129-143.PubMedPubMed CentralView ArticleGoogle Scholar
- Fang X, Wallqvist A, Reifman J: Modeling synergistic drug inhibition of Mycobacterium tuberculosis growth in murine macrophages. Mol Biosyst. 2011, 7: 2622-2636. 10.1039/c1mb05106g.PubMedView ArticleGoogle Scholar
- Becker SA, Feist AM, Mo ML, Hannum G, Palsson BO, Herrgard MJ: Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox. Nat Protoc. 2007, 2: 727-738. 10.1038/nprot.2007.99.PubMedView ArticleGoogle Scholar
- Pharkya P, Burgard AP, Maranas CD: OptStrain: a computational framework for redesign of microbial production systems. Genome Res. 2004, 14: 2367-2376. 10.1101/gr.2872004.PubMedPubMed CentralView ArticleGoogle Scholar
- Liu L, Agren R, Bordel S, Nielsen J: Use of genome-scale metabolic models for understanding microbial physiology. FEBS Lett. 2010, 584: 2556-2564. 10.1016/j.febslet.2010.04.052.PubMedView ArticleGoogle Scholar
- Schuster S, Dandekar T, Fell DA: Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering. Trends Biotechnol. 1999, 17: 53-60. 10.1016/S0167-7799(98)01290-6.PubMedView ArticleGoogle Scholar
- Becker SA, Palsson BO: Context-specific metabolic networks are consistent with experiments. PLOS Comput Biol. 2008, 4: e1000082-10.1371/journal.pcbi.1000082.PubMedPubMed CentralView ArticleGoogle Scholar
- Oberhardt MA, Goldberg JB, Hogardt M, Papin JA: Metabolic network analysis of Pseudomonas aeruginosa during chronic cystic fibrosis lung infection. J Bacteriol. 2010, 192: 5534-5548. 10.1128/JB.00900-10.PubMedPubMed CentralView ArticleGoogle Scholar
- Shlomi T, Cabili MN, Herrgard MJ, Palsson BO, Ruppin E: Network-based prediction of human tissue-specific metabolism. Nat Biotechnol. 2008, 26: 1003-1010. 10.1038/nbt.1487.PubMedView ArticleGoogle Scholar
- Fang X, Wallqvist A, Reifman J: Modeling phenotypic metabolic adaptations of Mycobacterium tuberculosis H37Rv under hypoxia. PLOS Comput Biol. 2012, 8: e1002688-10.1371/journal.pcbi.1002688.PubMedPubMed CentralView ArticleGoogle Scholar
- Narendra V, Lytkin NI, Aliferis CF, Statnikov A: A comprehensive assessment of methods for de-novo reverse-engineering of genome-scale regulatory networks. Genomics. 2010, 97: 7-18.PubMedPubMed CentralView ArticleGoogle 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-PubMedPubMed CentralView ArticleGoogle Scholar
- Covert MW, Xiao N, Chen TJ, Karr JR: Integrating metabolic, transcriptional regulatory and signal transduction models in Escherichia coli. Bioinformatics. 2008, 24: 2044-2050. 10.1093/bioinformatics/btn352.PubMedView ArticleGoogle Scholar
- Kim J, Reed JL: OptORF: Optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains. BMC Syst Biol. 2010, 4: 53-10.1186/1752-0509-4-53.PubMedPubMed CentralView ArticleGoogle Scholar
- Lee JM, Min Lee J, Gianchandani EP, Eddy JA, Papin JA: Dynamic analysis of integrated signaling, metabolic, and regulatory networks. PLOS Comput Biol. 2008, 4: e1000086-10.1371/journal.pcbi.1000086.PubMedView ArticleGoogle Scholar
- Kotte O, Heinemann M: A divide-and-conquer approach to analyze underdetermined biochemical models. Bioinformatics. 2009, 25: 519-525. 10.1093/bioinformatics/btp004.PubMedView ArticleGoogle Scholar
- Price ND, Reed JL, Palsson BO: Genome-scale models of microbial cells: evaluating the consequences of constraints. Nat Rev Microbiol. 2004, 2: 886-897. 10.1038/nrmicro1023.PubMedView ArticleGoogle 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: 1871-1878.PubMedView ArticleGoogle Scholar
- Liebermeister W, Klipp E: Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theor Biol Med Model. 2006, 3: 41-10.1186/1742-4682-3-41.PubMedPubMed CentralView ArticleGoogle Scholar
- Tran LM, Rizk ML, Liao JC: Ensemble modeling of metabolic networks. Biophys J. 2008, 95: 5606-5617. 10.1529/biophysj.108.135442.PubMedPubMed CentralView ArticleGoogle Scholar
- Jamshidi N, Palsson BO: Mass action stoichiometric simulation models: incorporating kinetics and regulation into stoichiometric models. Biophys J. 2010, 98: 175-185. 10.1016/j.bpj.2009.09.064.PubMedPubMed CentralView ArticleGoogle Scholar
- Smallbone K, Simeonidis E, Swainston N, Mendes P: Towards a genome-scale kinetic model of cellular metabolism. BMC Syst Biol. 2010, 4: 6-10.1186/1752-0509-4-6.PubMedPubMed CentralView ArticleGoogle Scholar
- Schwacke JH, Voit EO: Computation and analysis of time-dependent sensitivities in Generalized Mass Action systems. J Theor Biol. 2005, 236: 21-38. 10.1016/j.jtbi.2005.02.013.PubMedView ArticleGoogle Scholar
- Abbott DA, Knijnenburg TA, de-Poorter LM, Reinders MJ, Pronk JT, van-Maris AJ: Generic and specific transcriptional responses to different weak organic acids in anaerobic chemostat cultures of Saccharomyces cerevisiae. FEMS Yeast Res. 2007, 7: 819-833. 10.1111/j.1567-1364.2007.00242.x.PubMedView ArticleGoogle Scholar
- Moxley JF, Jewett MC, Antoniewicz MR, Villas-Boas SG, Alper H, Wheeler RT, Tong L, Hinnebusch AG, Ideker T, Nielsen J, Stephanopoulos G: Linking high-resolution metabolic flux phenotypes and transcriptional regulation in yeast modulated by the global regulator Gcn4p. Proc Natl Acad Sci USA. 2009, 106: 6477-6482. 10.1073/pnas.0811091106.PubMedPubMed CentralView ArticleGoogle Scholar
- Fournier ML, Paulson A, Pavelka N, Mosley AL, Gaudenz K, Bradford WD, Glynn E, Li H, Sardiu ME, Fleharty B, et al: Delayed correlation of mRNA and protein expression in rapamycin-treated cells and a role for Ggc1 in cellular sensitivity to rapamycin. Mol Cell Proteomics. 2010, 9: 271-284. 10.1074/mcp.M900415-MCP200.PubMedPubMed CentralView ArticleGoogle Scholar
- Lee MV, Topper SE, Hubler SL, Hose J, Wenger CD, Coon JJ, Gasch AP: A dynamic model of proteome changes reveals new roles for transcript alteration in yeast. Mol Syst Biol. 2011, 7: 514-PubMedPubMed CentralView ArticleGoogle Scholar
- Cleveland WS: Robust locally weighted regression and smoothing scatterplots. J Am Stat Assoc. 1979, 74: 829-836. 10.1080/01621459.1979.10481038.View ArticleGoogle Scholar
- Hinnebusch AG: Translational regulation of GCN4 and the general amino acid control of yeast. Annu Rev Microbiol. 2005, 59: 407-450. 10.1146/annurev.micro.59.031805.133833.PubMedView ArticleGoogle Scholar
- AbdulHameed MD, Chaudhury S, Singh N, Sun H, Wallqvist A, Tawa GJ: Exploring polypharmacology using a ROCS-based target fishing approach. J Chem Inf Model. 2012, 52: 492-505. 10.1021/ci2003544.PubMedView ArticleGoogle Scholar
- Piper P, Mahe Y, Thompson S, Pandjaitan R, Holyoak C, Egner R, Muhlbauer M, Coote P, Kuchler K: The pdr12 ABC transporter is required for the development of weak organic acid resistance in yeast. EMBO J. 1998, 17: 4257-4265. 10.1093/emboj/17.15.4257.PubMedPubMed CentralView ArticleGoogle Scholar
- Hatzimanikatis V, Bailey JE: MCA has more to say. J Theor Biol. 1996, 182: 233-242. 10.1006/jtbi.1996.0160.PubMedView ArticleGoogle Scholar
- Adiamah DA, Handl J, Schwartz JM: Streamlining the construction of large-scale dynamic models using generic kinetic equations. Bioinformatics. 2010, 26: 1324-1331. 10.1093/bioinformatics/btq136.PubMedView ArticleGoogle Scholar
- Bruck J, Liebermeister W, Klipp E: Exploring the effect of variable enzyme concentrations in a kinetic model of yeast glycolysis. Genome Inform. 2008, 20: 1-14.PubMedGoogle Scholar
- Teusink B, Passarge J, Reijenga CA, Esgalhado E, van der-Weijden CC, Schepper M, Walsh MC, Bakker BM, van-Dam K, Westerhoff HV, Snoep JL: Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. Eur J Biochem. 2000, 267: 5313-5329. 10.1046/j.1432-1327.2000.01527.x.PubMedView ArticleGoogle Scholar
- Kummel A, Panke S, Heinemann M: Putative regulatory sites unraveled by network-embedded thermodynamic analysis of metabolome data. Mol Syst Biol. 2006, 2: 2006.0034-PubMedPubMed CentralView ArticleGoogle Scholar
- Heinrich R, Rapoport TA: A linear steady-state treatment of enzymatic chains. General properties, control and effector strength. Eur J Biochem. 1974, 42: 89-95. 10.1111/j.1432-1033.1974.tb03318.x.PubMedView ArticleGoogle Scholar
- Kacser H, Burns JA: The control of flux. Symp Soc Exp Biol. 1973, 27: 65-104.PubMedGoogle Scholar
- Schwanhausser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, Chen W, Selbach M: Global quantification of mammalian gene expression control. Nature. 2011, 473: 337-342. 10.1038/nature10098.PubMedView ArticleGoogle Scholar
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.