Bioinformatics Advance Access originally published online on January 10, 2008
Bioinformatics 2008 24(5):704-710; doi:10.1093/bioinformatics/btn010
A feed-forward loop guarantees robust behavior in Escherichia coli carbohydrate uptake
Systems Biology Group, Max-Planck-Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Motivation: In Escherichia coli, the phosphoenolpyruvate: carbohydrate phosphotransferase system acts like a sensory element which is able to measure the flux through glycolysis. Since the output of the sensor, the phosphorylated form of protein EIIA, is connected to the activity of the global transcription factor Crp, the kinetic and structural properties of the system are important for the understanding of the overall cellular behavior.
Results: A family of mathematical models is presented, varying with respect to their degree of complexity (number of reactions that are taken into account, number of parameters) that show a structurally and quantitatively robust behavior. The models describe a set of experimental data that relates the output of the sensor to the specific growth rate. A central element that is responsible for the structural robustness is a feed-forward loop in the glycolysis, namely the activation of the pyruvate kinase reaction by a metabolite of the upper part of the glycolysis. The robustness is shown for variations of the measured data as well as for variations of the parameters.
Availability: MATLAB files for model simulations are available on http://www.mpi-magdeburg.mpg.de/people/kre/robust/ A short description of the files provided on this site can be found in the Supporting information.
Contact: kremling{at}mpi-magdeburg.mpg.de
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
In recent years, the set-up of mathematical models for cellular systems that describe metabolism, signal transduction and gene expression has become very popular and will lead to a better understanding of the underlying molecular processes. The knowledge on the detailed interactions between the components that are responsible for carbohydrate uptake in Escherichia coli is steadily increasing and current research on individual uptake systems like for glucose uptake via the phosphoenolpyruvate (PEP): carbohydrate phosphotransferase system (PTS) reveals new players that maybe play a role in control of these systems (Plumbridge, 1998). However, the knowledge on individual uptake systems is already rich and is used as a basis to set up mathematical models to describe and analyze the properties of the control circuits (e.g. see Bettenbrock et al., 2006; Mahadevan et al., 2002; Santillan and Mackey, 2004).
In previous reports (Bettenbrock et al., 2006, 2007; Kremling et al., 2007), we analyzed in detail a signal transduction pathway that senses the metabolic state of E.coli during carbohydrate uptake and processes the signal to activate Crp. Crp is a global transcription factor involved in the expression of a large number of genes, responsible for carbohydrate uptake and chemotaxis. A key element in this process is the PTS shown in Figure 1. The PTS is a transport and a sensory system at the same time. In a sequence of four reactions, a phosphoryl group is transferred from metabolite PEP to protein EIIACrr (EIIA is used further in the text), the output of the sensory system. For example, in case of glucose, the phosphoryl group is afterwards transferred the the actual transport protein EIICBGlc and then to the incoming sugar. Interestingly, a relationship between the specific growth rate µ and degree of phosphorylation of EIIA could be seen in various growth situations of the wild type strain growing on single substrates like glucose, lactose or glycerol, and also for growth on mixtures of substrates (see Bettenbrock et al., 2007 and Hogema et al., 1998).
|
Often sensor elements can be regarded as logic elements that process external stimuli into intracellular signals. As an example, a NOT element with high input will result in a low response of the output. Circuits representing different logic elements are mainly found in signaling cascades of higher cells. Here, we present a logic element that can be found in bacterial metabolism: high fluxes through the glycolysis, corresponding to high growth rates result in a low degree of phosphorylation of EIIA. This is surprising, since, assuming a linear reaction chain, high fluxes result in high concentrations of the metabolites in the pathways, based on the (normally) monotone dependency of the reaction rate from the substrate concentrations. The PTS together with the glycolysis can now be seen as an element that allows the transformation of high fluxes into a low-metabolite concentrations.
Robustness is the insensitivity of a selected characteristic (time course of a component, steady-state characteristics, network function to sustain growth (Stelling et al., 2002), adaption precision (Barkai and Leibler, 1997) with respect to changes of external or internal perturbations (different environmental conditions, mutations, altered kinetic parameters or altered model structures). For the contribution at hand, we define a set of structural and quantitative robust mathematical models as models that fulfill the following conditions:
- the models are quantitative, that is, they describe experimental data (time course data or steady-state characteristics) representing a cellular function with a given accuracy;
- the models of the family differ in the number of components, the number of reactions, the number of regulatory pattern and/or in the choice of the kinetic expressions for the reaction rates.
|
2 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Figure 2 shows a scheme of the biochemical network that is responsible for metabolism of carbohydrates. In general, substrates enter glycolysis at different nodes. The scheme in Figure 2, left side, considers substrates that feed into glucose 6-phosphate, the first glycolytic metabolite. The scheme covers the central reactions of carbohydrate metabolism. The state variables are Glc6P (glucose 6-phosphate), TP (triose phosphate), PEP (phosphoenolpyruvate), Prv (pyruvate), and EIIA / EIIAP. EIIA / EIIAP represent all the PTS proteins. Since the reactions of the PTS are very fast in comparison to glycolytic reactions or gene expression (Kremling et al., 2004), the individual reactions of the PTS can be lumped together and can be described by a single equilibrium constant Kpts. Rates rup_npts and rup_pts represent uptake of either a non-PTS sugar or of a PTS sugar, respectively. Rates rppp/bio, rppp, rtca/bio represent fluxes from glucose 6-phosphate into pentose phosphate pathway and biosynthesis, flux from pentose phosphate pathway back to glycolysis and drain to TCA and biosynthesis from PEP, respectively. Fluxes through glycolysis are represented by rgly, reno and rpyk and it is assumed that the respective enzyme concentrations depend on the growth rate µ. Rate rpts represents the rate through the PTS. The equations for the state variables are summarized in the Supporting information.
|
A feed-forward loop, the activation of the pyruvate kinase by a metabolite from the upper part of the glycolysis is described in the literature (Waygood and Sanwal, 1974) (right side of Fig. 2). The activator of the pyruvate kinase is fructose 1,6-bis-phosphate. In the models introduced below, fructose 1,6-bis-phosphate is not included as a state variable. Therefore, it is replaced by glucose 6-phosphate. This is justified since in the upper part of the glycolysis the drain to anabolism is very small and it can be expected that the steady-state values do not differ very much between glucose 6-phosphate and fructose 1,6-bis-phosphate.
2.1 Sensory system
At first, the sensory system is considered and a relationship for the output of the PTS, phosphorylated EIIA, has to be derived. If the flux distribution at PEP is considered for the case that the PTS is present but not involved in uptake (e.g. growth on non-PTS sugars like glucose 6-phosphate, glycerol, etc.) the reaction rate of the reversible reaction rpts has to be zero:
|
| (1) |
. Taking into account that the PTS proteins are either phosphorylated or not, the overall concentration EIIA0 is introduced:
|
| (2) |
|
| (3) |
. In case of an active PTS, i.e. the carbohydrate is taken up and phosphorylated by the PTS, the general structure of Equation 3 is also valid. The following steady-state equation will hold for PTS substrates:
|
| (4) |
|
| (5) |
|
The PEP and pyruvate concentrations in a cell are difficult to measure. The techniques established for the measurement of metabolites generally generate data with high errors and reliable data about the variation of PEP and pyruvate concentrations with increasing growth rate are lacking. The PEP to pyruvate ratio in a cell can be measured indirectly via the phosphorylation level of EIIA and this ratio has been shown to decrease with increasing growth rate. The networks shown in Figure 2 lead to two different scenarios that are depicted in Figure 3: high growth rates require high fluxes and lead to increased fluxes through the pyruvate kinase. In a network without further control both, the PEP and the pyruvate concentrations have to increase with increasing growth rate (see Fig. 3) to match this requirement. In this case, the PEP to pyruvate ratio is very sensitive to small fluctuations in the metabolite concentrations and is most probably difficult to control. On the contrary, in a network that is controlled via the feed-forward loop shown in Figure 2 (right plot) high fluxes can be realized by lowering the PEP concentration. In this case, the pyruvate concentration increases with increasing growth rate while the PEP concentration decreases. The high flux through the pyruvate kinase, in this case, is realized by the activation by glucose 6-phosphate. This scenario is much less fragile to small fluctuations or uncertainties in the metabolite concentrations, since the ratio of a decreasing metabolite (PEP) and an increasing metabolite (pyruvate) always decreases.
In the following, different model variants are introduced that correspond to two cases: uncontrolled network or controlled network with feed-forward loop.
2.2 Model variants
The model variants that are used for this study range form detailed to simple, but are all based on the available biological knowledge. The models differ (i) in the kinetic expression for rgly, rpyk and rpdh, (ii) in the dependency of the enzyme concentration on the growth rate and (iii) in the incorporation of the drain into biosynthesis. With this approach, models of different complexity are generated and analyzed with respect to model verification by experimental data.
Model 1 considers all dependencies shown in Figure 2 (left plot). Fluxes rppp/bio, rppp and rtca/bio are 40, 20 and 25% of the uptake rate (Holms, 1996). Levels of enzymes depend on the activity of the transcriptional/translational machinery and on the activity of transcription factors. This results in different concentrations of enzymes if the whole range of growth rates is considered. For example, Seeto et al. (2004) report on the dependency of PtsG from the diluation rate during continuous cultivation. Experimental data for catalytic enzymes like in the glycolysis are, however, not available. To take into account possible dependencies, a simple relationship was included here: the dependency of scaled glycolytic enzyme concentrations e/e0 from growth rate µ is as follows:
|
| (6) |
The models also differ in the choice of the kinetics for the glycolysis reaction rgly, rpyk and rpdh. Here, mass action law or Michaelis–Menten kinetics are used. Model 1 is the uncontrolled network. The variants of Model 1 that are used i the analysis are summarized in Table 1.
|
Model 2 represents the controlled network. In the controlled network, the pyruvate kinase reaction is controlled by a feed-forward loop. Different rate laws are used to describe the reaction rates; e.g. one choice for pyruvate kinase is based on recent publications that use a Monod–Wyman–Changeux kinetics (Bettenbrock et al., 2006; Chassagnole et al., 2002). Variants of Model 2 that are used for analysis are summarized in Table 1.
2.3 Parameter estimation and model assessment
The objective function to fit the parameters was formulated as a ordinary least square problem and a standard gradient-based algorithm as it is provided by MATLAB was used for solving. In the Section 1, we suggest to use a measure that allows to assess quantitatively the results of the parameter estimation. Since the measurement error of the measurements [experimental data used are described in detail in (Bettenbrock et al., 2006, 2007)] can hardly be determined we calculate the estimated standard deviation [or residual mean square (Montgomery et al., 2001)] 
. Finally, we relate 
on the overall concentration EIIA0 to get a % value:
|
| (7) |
i, N the number of data points and n the number of parameters that were used in the estimation. Experiments were performed with different substrates and substrate combinations under different conditions (Bettenbrock et al., 2006, 2007) and N = 45 data points are available. The number of estimated parameters is n = 4, i.e. kgly, kpyk, kpdh and kpts are estimated.
The model with the best fit is Model 2b, that is, the model with the highest complexity (Monod–Wyman–Changeux kinetics for the pyruvate kinase, drain fluxes to monomers as well as growth dependent enzyme concentrations are considered). In general, all variants of Model 2 have better results than the models of class Model 1, except for Model 1c. Interestingly, Model 2f with the simplest structure and the simplest rate laws but taking into account the feed-forward loop also reaches a good 
value. Figures 5 and 6 compare some of the model variants with experimental data.
2.4 Model analysis
In this section, the sensitivity of the objective function is analyzed with respect to (i) variations of the measured data and (ii) with respect to model parameters. For the analysis one representative of the model class M1 (uncontrolled network) and one representative of model class M2 (controlled network) were considered. The models are comparable since their structure is the same except the feed-forward loop. It is expected that other model variants behave similarly than the representatives.
2.4.1 Variations of the measured data
Recently, we applied a statistical procedure—the bootstrap method (DiCiccio and Efron, 1996; Efron and Tibshirani, 1993)—to determine parameter uncertainties for nonlinear systems (Joshi et al., 2006). The method surmounts the theoretical limitations (e.g. the Fisher-Information matrix gives only a lower bound for the parameter variances in case of systems that are nonlinear in parameters) by assessing the uncertainties in statistics with data from finite samples. Like a Monte-Carlo method, the bootstrap method uses stochastic elements and repeated simulations to analyze the properties of the system under consideration.
Briefly, the analysis is performed in such a way, that an initial set of experimental data S is used as a data base. Performing parameter estimation result in a first set of parameters and value to assess the model quality. Due to measurement errors the repetition of the experiment leads to a slightly different set of data S1 and therefore to a different set of estimated parameters and value. The bootstrap approach now uses a large set of B-times replicated experimental data S1, S2, S3, ... ,SB to calculate statistical properties of the resulting distribution of the (re)-estimated set of parameters and values. Formal we look for
|
| (8) |
with respect to modification of the measured data 
. The values for CR95 can directly be read off from the distributions shown in the following figures. Figure 7 shows the results of the bootstrapping with respect to measurement variations for model M1 and M2. For each of the 2000 runs, every single data point was modified by adding a random noise (normal distribution) of 2%. As can be seen, model variant M2d representing the controlled network shows a narrower distribution of the values.
2.4.2 Variations of parameters
An analogous approach is applied to determine the influence of parameter variations p on the values:
|
| (9) |
], the model parameters that were used for parameter estimation (kgly, kpdh, kpts and kpyk) are altered by adding a random number again from a normal distribution (10%), simulating the system, and calculating with Equation 7. Figure 8 reveals that the model representing the controlled system shows a narrower distribution. |
3 CONCLUSIONS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Robustness in cellular systems has been discussed very frequently in the last years. A prominent example is the bacterial chemotaxis where an integral feedback loop is responsible for precise adaptation of the system with respect to internal perturbations.
The signal flow, responsible for protein synthesis in the E.coli carbohydrate uptake network has been under investigation for a long time and several players have been described in the literature. Besides local control by carbohydrate specific regulators, the global regulator Crp is involved in transcription initiation for almost all genes in the network. Experimental data revealed a relationship between the specific growth rate of E.coli and the output of the sensor phosphorylated EIIA of the signaling pathway. The sensor measures the flux through the central pathways. This is realized by a network structure that transforms a high flux into a low response and can be seen therefore as a logic element with NOT function. A number of model variants are thinkable that allow a quantitative description of the experimental data. Therefore, a number of different models were set up and analyzed.
Model development is always a competition between a realistic description, that is based on the available knowledge, and a reduced or simplified description, that takes into account only the most important characteristics. In this contribution, the models presented show a different degree of complexity based on the kinetic expression for the single rates or the number of reactions that are taken into account. The idea behind this is to show that available experimental data can be described with good accuracy not only by a single model structure but by a whole class of models. Here, a number of different model structures, taking into account different flux distributions, kinetic expressions and possible dependencies of the enzyme concentrations on the growth rate are set up and investigated. The analysis reveals that only those model variants that include a special motif, a feed-forward loop, can describe the data with high accuracy. This feature is named quantitative robustness, since the reproduction of experimental data—here a characteristic curve—is required. The minimal value is achieved with model variant M2b ( = 13.19%) and the maximum with model variant M2f ( = 15.2%), indicating that the model variants show nearly equal accuracy. For the models without feed-forward loop the difference is much larger (the values are between 14.9% and 23.6%, see also Fig. 4).
|
|
non-PTS carbohydrates and 
PTS carbohydrates from Bettenbrock et al. (
non-PTS carbohydrates and 
PTS carbohydrates from Bettenbrock et al. (
|
|
|
The analysis of the circuit reveals that the feed-forward loop is a robust element. Small variations or disturbances will affect the function of the PEP/pyruvate ratio only marginally (Fig. 3). We expected that with model variants M2 it is possible to describe also slightly different experimental data with higher accuracy than with model variants M1. Therefore, a bootstrap approach was performed and the analysis reveals that indeed model variants M2 give better results than model variants M1 (Fig. 7). The 95% interval for model variant M1 CR95 = 16% while for model variant M2 CR95 = 10.28%. To assess the influence of the parameters, the four parameters that were estimated are randomly modified and
values are calculated. For Model M1d the 95% interval CR95 = 9.36% while for model M2d CR95 = 7.17%. From the results it could be concluded that model variants M2, including the feed-forward loop are more appropriate to describe the available data because they show better structural and quantitative characteristics than models without the loop. |
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
A.K. and K.B. are funded by the FORSYS initiative from the German Federal Ministry of Education and Research (BMBF).
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Thomas Lengauer
Received on July 5, 2007; revised on December 10, 2007; accepted on January 6, 2008
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 RESULTS3 CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Barkai N, Leibler S. Robustness in simple biochemical networks. Nature (1997) 387.
Bettenbrock K, et al. A quantitative approach to catabolite repression in Escherichia coli. J. Biol. Chem (2006) 281:2578–2584.
Bettenbrock K, et al. Analysis of the correlation between growth rate, EIIACrr (EIIAGlc) phosphorylation levels and intracellular cAMP levels in E.coli K-12. J. Bacteriol (2007) 189:6891–6900.
Chassagnole C, et al. Dynamic modeling of the central carbon metabolism of E.coli. Biotech. Bioeng (2002) 79:53–73.[CrossRef]
DiCiccio TJ, Efron B. Bootstrap confidence intervals. Stati. Sci (1996) 11.
Efron B, Tibshirani. RJ. An Introduction to the Bootstrap. (1993) Boca Raton, Florida: Chapman and Hall.
Hogema BM, et al. Inducer exclusion in Escherichia coli by non-PTS substrates: the role of the PEP to pyruvate ratio in determining the phosphorylation state of enzyme IIAGlc. Mol. Microbiol (1998) 30:487–498.[CrossRef][Web of Science][Medline]
Holms H. Flux analysis and control of the central metabolic pathways in Escherichia coli. FEMS Microbiol Rev (1996) 19:85–116.[CrossRef][Web of Science][Medline]
Joshi M, et al. Exploiting the bootstrap method for quantifying parameter confidence intervals in dynamical systems. Metab. Eng (2006) 8:447–455.[CrossRef][Medline]
Kremling A, et al. Time hierarchies in the Escherichia coli carbohydrate uptake and metabolism. BioSystems (2004) 73:57–71.[CrossRef][Medline]
Kremling A, et al. Analysis of global control of Escherichia coli carbohydrate uptake. In: BMC Syst Biol. (2007) 1:42.
Mahadevan R, et al. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys. J (2002) 83:1331–1340.[Medline]
Mangan S, et al. The coherent feedforward loop serves as a sign-sensitive delay element in transcription networks. J. Mol. Biol (2003) 334:197–204.[CrossRef][Web of Science][Medline]
Montgomery DC, et al. Engineering Statistics. (2001) New York: John Wiley and Sons.
Plumbridge J. Expression of ptsG, the gene for the major glucose PTS transporter in E.coli, is repressed by Mlc and induced by growth on glucose. Mol. Microbiol (1998) 29:1053–1063.[CrossRef][Web of Science][Medline]
Santillan M, Mackey. MC. Influence of catabolite repression and inducer exclusion on the bistable behavior of the lac operon. Biophys. J (2004) 86:1282–1292.[Web of Science][Medline]
Seeto S, et al. The multifactorial influences of RpoS, Mlc and cAMP on ptsG expression under glucose limited and anaerobic conditions. Res. Microbiol (2004) 155:211–215.[Medline]
Stelling J, et al. Metabolic network structure determines key aspects of functionality and regulation. Nature (2002) 420:190–193.[CrossRef][Medline]
Waygood EB, Sanwal. BD. The control of pyruvate kinase of Escherichia coli. J. Biol. Chem (1974) 249:265–274.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  O. Kotte and M. Heinemann A divide-and-conquer approach to analyze underdetermined biochemical models Bioinformatics, February 15, 2009; 25(4): 519 - 525. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||