Sepideh
Dolatshahi
,
Luis L.
Fonseca
and
Eberhard O.
Voit
*
Department of Biomedical Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA. E-mail: sepideh@gatech.edu; llfonseca@gatech.edu; eberhard.voit@bme.gatech.edu
First published on 6th November 2015
This article and the companion paper use computational systems modeling to decipher the complex coordination of regulatory signals controlling the glycolytic pathway in the dairy bacterium Lactococcus lactis. In this first article, the development of a comprehensive kinetic dynamic model is described. The model is based on in vivo NMR data that consist of concentration trends in key glycolytic metabolites and cofactors. The model structure and parameter values are identified with a customized optimization strategy that uses as its core the method of dynamic flux estimation. For the first time, a dynamic model with a single parameter set fits all available glycolytic time course data under anaerobic operation. The model captures observations that had not been addressed so far and suggests the existence of regulatory effects that had been observed in other species, but not in L. lactis. The companion paper uses this model to analyze details of the dynamic control of glycolysis under aerobic and anaerobic conditions.
This article and its companion^{1} describe this process of connecting snapshots into cohesive storylines within the context of a complex task that the bacterium Lactococcus lactis must solve. Namely, the bacterium has to survive periods of starvation and at the same time be optimally positioned to take up substrate once it becomes available again, which is not a trivial matter, as we will discuss. Expressed differently, it must reach a “ready-to-respond” state while substrate is running out. The goal of these two papers is to decipher through dynamical modeling how the cells coordinate a complex control program to ensure this state. The first paper focuses on the construction of the model, while the second paper^{1} utilizes this model to shed light on the complex regulatory program, with which the organism controls glycolysis and ensures a viable ready-to-respond state.
The construction of a dynamic model of the type envisioned here has become possible with the recent availability of comprehensive time series data that characterize the dynamics of genomic, proteomic, metabolic, and physiological responses. Of particular interest for biochemical pathway modeling is the availability of relatively dense time profiles of metabolites through measurement techniques such as mass spectrometry (MS^{2,3}), nuclear magnetic resonance (NMR^{4}), protein kinase phosphorylation, or mass cytometry (CyTOF^{5}). These data contain valuable, but implicit, information about the structure and dynamics of the biological system under study and permit the use of top-down modeling approaches. These approaches consist of minimizing the residual error between the measured data, i.e., the time profiles, and the assumed model, which typically consists of a system of nonlinear ordinary differential equations (ODE) that are to be parameterized.^{6}
In an effort to build a kinetic-dynamic model using time-series data, the most challenging steps are the identification of suitable mathematical formats for all flux representations and the estimation of their unknown kinetic parameters. Central to these tasks is the methodology of dynamic flux estimation (DFE).^{7}
In contrast to previous models, which captured the dynamics of single datasets (e.g., ref. 8 and 9), this paper describes a single aggregate model that combines three datasets for different input glucose concentrations. This single model enables the assessment of several intriguing observations, which to date had not been explained convincingly.
A particular observation that has been puzzling for a long time is the following: if glucose is supplied to the medium in different amounts, one would expect a corresponding increase in the peak level of the downstream metabolite fructose-1,6-bisphosphate (FBP). However, the data clearly show that the peak concentration is largely independent of the external glucose concentrations. Specifically, comparing the different series of experimental results with increasing glucose concentrations, the FBP accumulation shows a progressively more noticeable plateau, whose duration, but not height, varies with substrate availability (Fig. 1). A previous model for aerobic conditions^{9} was able to fit glycolytic data for a single glucose concentration of 20 mM, but extrapolating the model toward different glucose inputs, such as 40 or 80 mM, predicted almost a doubling or quadrupling of the FBP peak height, which is in obvious contrast to the experimental measurements under anaerobic conditions. Here, we analyze the situation in detail and extract constraints that the model and its parameters must satisfy to yield a good fit to the data.
Fig. 1 FBP data for three external glucose concentrations under anaerobic conditions. The peak level appears to be independent of the available substrate concentrations. One also notes that FBP does not vanish completely and instead maintains some residual concentration. Data summarized from ref. 12 and 13. |
Throughout the analysis here and in the companion paper, such sometimes non-intuitive observations will be discussed, and the roles of numerous precisely timed regulatory mechanisms will be explained through controlled simulations with and without such mechanisms.
The measured time series include the external concentrations of glucose and of the end product lactate in the medium, along with several of the more abundant intermediate metabolites, namely glucose 6-phosphate (G6P), fructose 1,6-bisphosphate (FBP), phosphoenol pyruvate (PEP), 3-phosphoglycerate (3PGA). The time series reflect three experiments performed with 20, 40, or 80 mM of glucose input, respectively. Additional time series are available for some ubiquitous cofactors, such as NAD^{+} and NADH, which are detectable with somewhat limited detection capability due to the required NMR acquisition time.^{12} Finally, the level of NTP was measured using ^{31}P-NMR. The datasets used for modeling are summarized in Table 1.
Experiment | Condition | Technique | Metabolites measured | Time interval |
---|---|---|---|---|
(1) 20 mM [6-^{13}C] glucose |
Anaerobic
pH = 6.5 |
^{13}C-NMR | glc, G6P, FBP, 3PGA, lac | 30 s |
(2) 40 mM [1-^{13}C] glucose and [5-^{13}C] nicotinate |
Anaerobic
pH = 6.5 |
^{13}C-NMR
^{31}P-NMR |
glc, FBP, 3PGA, PEP, lac, NAD^{+}, NADH
ATP, P_{i} |
2.2 min
2.75 min |
(3) 80 mM [1-^{13}C] glucose and [5-^{13}C] nicotinate |
Anaerobic
pH = 6.5 |
^{13}C-NMR
^{31}P-NMR |
glc, FBP, 3PGA, PEP, lac, NAD^{+}, NADH
ATP, P_{i} |
2.2 min
2.75 min |
In Experiment 1, 20 mM of [6-^{13}C] glucose were supplied to the cell suspension, and time series of concentrations were recorded for glucose, G6P, FBP, 3PGA, and lactate with a resolution of 30 seconds. In Experiments 2 and 3, cells were supplied beforehand with labeled [5-^{13}C] nicotinic acid, a precursor of NADH, thus ensuring this pool to be 100% labeled.^{12} To start the in vivo NMR experiments, cells were supplied with 40 and 80 mM of [1-^{13}C] glucose, respectively, and time series of glucose, FBP, PEP, 3PGA, lactate, NAD^{+} and NADH were recorded. The time resolution in these experiments was 2.2 minutes, which was needed to accommodate the determination of NAD^{+} and NADH with a reasonable signal-to-noise ratio. In a separate, comparable experiment, cellular metabolism was investigated with ^{31}P-NMR, which allowed the measurement of NTP (mostly ATP), P_{i} and pH with a time resolution of 2.75 min. Usage of [1-^{13}C] glucose prevents the determination of G6P in these datasets, due to the similarity in chemical shifts of [1-^{13}C] glucose and [1-^{13}C] G6P. However, some measurements of G6P are available from experiments with 20 mM [6-^{13}C] glucose. Although the data described above exhibit experimental noise and are incomplete, with some metabolites or time points missing in each set, they are as good as a modeler can presently hope for. Duplicates are available for Experiments 1 and 2. These replicate datasets have a highly similar behavior to the modeled datasets, which shows that the data obtained with these NMR methods are highly reproducible. The additional datasets are not used for the modeling for the same reason of similarity and the fact that they contain time series data for fewer intermediate metabolites.
Under the given experimental conditions, cells are not able to synthesize ATP (nor ADP) de novo, because they are suspended in phosphate buffer and supplied only with glucose, which is insufficient for L. lactis to grow. Once the culture is harvested and washed prior to the NMR experiment, the total amount of ATP + ADP remains constant. This fact enables us to infer ADP from the ATP data where these are available (for Experiments 2 and 3), which is beneficial for the initial parameter estimation. By the same token, the cells do not synthesize proteins and they do not alter the expression of their genes.
In all datasets, the concentrations of 3PGA and PEP are highly covariant, due to the fact that these two metabolites can be converted into 2PGA by two enzymatic steps (phosphoglycerate mutase and enolase). These reactions are fast and reversible with equilibria that favor 3PGA and PEP,^{14–16} thus maintaining 2PGA in a concentration range below the in vivo^{13}C-NMR detection level. Since this covariance is maintained even during high glycolytic flux, we decided to aggregate these two pools into one dependent variable (X_{4}). Nonetheless, we are still able to calculate the concentration of each intermediate from the constant proportionality of ∼0.6444 (3PGA/PEP).
The data seem to indicate an apparent absence of PEP and FBP at the beginning of the experiment. However, these metabolites are present, but just not labeled initially and therefore missed by the NMR detection. The concentrations of these metabolites were measured in a control experiment,^{11} supporting the reasonable assumption that the cells re-enter a state of starvation at the end of the experiment and that the residual values of now labeled PEP and FBP constitute a state that is similar to the state at the beginning of the experiment. This conjecture is further supported by a tandem experiment in which glucose was presented twice after periods of starvation.^{10}
The structure of the pathway at hand is presented in Fig. 2. A dynamic model for a similar type of experiments with 20 mM external glucose, but under aerobic conditions, was developed previously.^{9} We started out with this model attempting to adjust it toward anaerobic conditions by including additional regulatory mechanisms from the biochemical literature, accounting for ubiquitous metabolites including NAD^{+}, NADH, ATP, and ADP as dependent variables, and making other necessary adjustments to the pathway based on biological literature. Unfortunately, the final task of parameter estimation did not succeed when the three datasets were to be captured simultaneously. Thorough analysis suggested a revisiting of the assumptions that were made originally regarding the pathway regulation and of the subtle, but genuine differences between operation under aerobic and anaerobic conditions. The amendments and alterations that became necessary will be explained in the following sections.
We use for this identification step dynamic flux estimation (DFE).^{7} In a nutshell, DFE works in the following manner. First, the time courses of all metabolites are smoothed, so that slopes can be estimated. These slopes are substituted for the derivatives on the left-hand sides of the ordinary differential equations (ODEs) at many time points, thereby replacing each ODE with a set of algebraic equations. At each time point, the algebraic equations can be solved by matrix inversion, with the help of pseudo-inverses, or with auxiliary methods.^{21–24} Combining these solutions yields graphical or numerical representations of all fluxes, either plotted against time or against the variables on which they depend. These representations can secondarily be fitted with functions that might be appropriate, which directly permits an assessment of the appropriateness and quality of the chosen functions.
The main attraction of DFE is the fact that this method does not presuppose a functional form for any of the flux representations. This feature allows us to test in an objective manner whether particular functions, such as power-laws, Michaelis–Menten rate laws, or Hill functions, are capable of appropriately modeling a specific flux, or if other formulations should be considered. Importantly, careful analyses of all fluxes in this manner may suggest the existence of regulatory signals that had been missing from the assumed pathway structure. Such a suggestion corresponds to a novel hypothesis that is in principle testable with lab experiments and may lead to biological discoveries. An example is the glucose uptake process, which the companion paper^{1} discusses in detail. In this case, slopes can directly be obtained through numerical differentiation of the uptake profile. The result of the first phase of DFE, which is independent of any model assumptions, reveals that this process must be regulated and suggests candidate regulators. These insights are used to refine the model construction.
In addition to its diagnostic capacities, DFE allows for a much more efficient parameter estimation strategy in terms of computation costs that are associated with the integration of ODEs and with global optimization. Namely, the parameters can be estimated one flux at a time. This decoupling is very advantageous, as one uses explicit functions, thereby avoiding—or at least reducing—the need to integrate ODEs numerous times, and because there is no risk of error compensation among terms or equations.
The first, model-free phase of DFE requires the estimation of slopes from the metabolic time course data. Several methods have been proposed for this purpose (e.g., see discussion in ref. 25); as a further alternative, a smoothing algorithm was introduced in ref. 26 based on the wavelet transforms that ensures conservation of mass.
Not all processes in the model can be directly assessed with DFE because; for instance, some time courses are missing. Nonetheless, DFE provides a beneficial initial step that significantly constrains the subsequent parameter estimation.
Our generic strategy thus is the following. As far as feasible, we use DFE to identify candidates for possibly missing regulatory effects. Next, we computationally select reasonable functional forms in a process that is ideally guided by biology known from the literature. Then we use DFE to obtain starting points for as many parameters as possible. For the remaining parameters, we use multi-start and randomized search algorithms, implemented with generous search intervals, to make sure that the candidate parameter sets are not unduly restricted. We do not impose limits on the rate constants, except that they have to be non-negative, and allow the kinetic orders to vary within boundaries of [−3,3], which are considered in the field as rather wide. Once a working parameter set is obtained, we use evolutionary and randomized optimization methods, followed by local searches, to attain the best parameter set according to the criterion of minimizing the sum of squared errors between data and model simulations. Finally, a search in the vicinity of the optimal solution will indicate the sloppiness of the solution, and large-scale random searches outside the vicinity indicate whether the solution is more or less unique.
One should note that the model construction process is not always cut and dry. In particular, some of the insights obtained in the companion paper fed back into the estimation process, as they constrained the parameter space. More generally, the state of the art is at this point not advanced enough to permit a fully automated parameter estimation of complex models.
In our specific case, several factors render the estimation of system parameters particularly challenging. First, the regulatory structure and reaction mechanisms of the model are a priori unknown, and their identification is by itself non-trivial, even if aided by DFE. Second, the error surface is embedded within a large parameter space of 44 dimensions. This surface is complicated, and one has to expect numerous local minima that severely confound steepest descent, evolutionary, or randomized optimization techniques that are typically used for parameter estimation. Third, datasets for some of the intermediate metabolites are not available, due to experimental limitations and, in particular, of the NMR-technique, which include higher detection limits for the various metabolites than would be desirable. Also, the ^{13}C- and ^{31}P-NMR experiments were executed separately and therefore both exclude some of the dependent variables. Missing data include time series for pyruvate for all three datasets, G6P for Datasets 2 and 3, and ATP, NAD^{+}, and NADH for Dataset 1, as well as F6P, G3P, DHAP, 1,3GBP, 2PGA, and ADP for all datasets. Fourth, stiffness and other numerical issues associated with the model equations can lead to substantial algorithmic difficulties while integrating the differential equations. While these issues may not even be present in the ultimate, optimized model, they are frequently encountered when an estimation algorithm determines inopportune combinations of settings during its scanning of a high-dimensional parameter space.^{27} Finally, even if a good fit has been determined, one cannot be sure of its uniqueness and must assume that the system might have a certain degree of sloppiness.^{28–32} In the present case, incorporating three datasets in conjunction with DFE can alleviate the identifiability issues^{30} and help reveal regulatory information about the pathway.
A generic reaction v_{ik} within the ith equation of a BST model is represented as , where the rate constant α_{ik} and the kinetic orders h_{ikj} are fundamental characteristics. The rate constant describes the turn-over of the process, while the kinetic orders quantify the strengths with which reactants and regulators affect the process. These parameters need to be estimated to provide a full description of the metabolic pathway under study.
(1a) |
(1b) |
(1c) |
(1d) |
(1e) |
(1f) |
(1g) |
(1h) |
(1i) |
(1j) |
(1k) |
(1l) |
(2a) |
(2b) |
(2c) |
(2d) |
(2e) |
(2f) |
(2g) |
(2h) |
(2i) |
Several options are available to address the initial brief rise in PTS. First, one could attempt to identify the true mechanisms leading to the initial increase in glucose uptake, which corresponds to the sigmoidal shape of the glucose concentration curve. For instance, one could explain this observation with the fact that the experimental set-up causes slight delays, which could affect the uptake profile. It has also been argued^{9} that the cells might recover from starvation with some variability, which could be due to cell-to-cell variation in the speed of glucose uptake or slight differences in glucose availability to individual cells, or could be caused by the process of mixing of glucose throughout the medium or other extraneous factors. Indeed, it was shown with simulations that a narrowly distributed uptake profile can directly convert the monotonic trend in glucose consumption, as predicted by the model structure, into a sigmoidal trend, as it is observed.^{9} As an alternative, one could try to identify some fitting function without being constrained by mechanistic considerations.^{41} One could furthermore use the data directly as (“off-line”) inputs instead of representing them functionally in the model.^{42} Finally, one could imagine that insufficient amounts of FBP accumulate during the first two minutes, thereby limiting the material flux into the pool of PEP and 3PGA. As a consequence, the initial PEP concentration might not be able to produce the high level of v_{1} that we directly computed through numerical differentiation of the glucose consumption profile. This short-term discrepancy can be resolved when we account for a small auxiliary permease flux (v_{P}), which has been observed in this strain of L. lactis,^{43} but was never used in earlier models.
Given the variety of speculations, we decided on a hybrid option, where we use a black-box adjustment function for the first two minutes and dynamically model glucose uptake afterwards, starting at t = 2 [min]. Specifically, we multiply a term of the form to the modeled PTS flux. This term, which is borrowed from physics and engineering,^{44} where it typically represents a delayed response of a linear time-invariant system to a constant input starting at time zero, does not convey biological meaning. It is simply useful since it exponentially approaches 1 and loses its effect after a short period of time. For instance, at 3T minutes, it is equal to 1 − e^{−3} = 0.9502. Thus, the term is ineffective for most of the experimental period. Different time constants T were allowed for different experiments.
The analysis led to several slight amendments of the original diagram. First, we observed that the actually measured total carbon mass decreases over time for all three datasets, with about 7–10% of the mass being unaccounted (see Fig. 6A in a later section), in spite of the fact that, outside glycolysis, the organism is metabolically inactive under the given conditions.^{10} The loss is biologically not very significant, but it is immediately inconsistent with the model structure of the initial diagram, where glucose is the only input substrate and lactate is the only output. The most reasonable option for remedying this inconsistency is the addition of a minor efflux out of one or more metabolite pools. A good candidate is pyruvate, because several reactions may use pyruvate as a substrate and thereby be responsible for the diversion of material.^{45} Because some of these effluxes potentially affect the NAD^{+}/NADH balance, we included two types of leakage from pyruvate, one with and one without the consumption of NADH. Optimization determined the very low capacities of these fluxes. We also considered alternative locations of leakage, such as PEP and G6P, but did not find them beneficial under the given experimental conditions. Other amendments are discussed in the following sections.
The model contains twelve dependent variables. Six represent the main metabolites glucose, G6P, FBP, the aggregated pool of 3PGA and PEP, pyruvate, and lactate, four variables represent the cofactors ATP, ADP, NAD^{+}, and NADH, and the remaining two represent temporary buffers.^{46} The buffers were used for the following reason. When we studied both carbon flow and electron flow, we observed a temporal mismatch that required some NADH to be consumed but later returned to the pathway. We rationalized this observation with the possibility that one or more pools must exist in connection with some glycolytic intermediates that allow NADH to be consumed and later returned. These reactions could include processes such as the mannitol-1P dehydrogenase or glycerol-3P dehydrogenase steps. The fact that there are several possible pools might explain why no other metabolite was actually measured by the NMR, because if the missing mass is distributed among several pools, none would have a concentration high enough to be detectable by the NMR. We therefore allowed for two buffers to exist in the model, one out of G6P and one out of pyruvate. No parameter set was ever obtained, by optimization, which used the buffer out of G6P, which suggested removal of this buffer. The best data consistency was always obtained with the two buffers presented, one NADH dependent and one NADH independent. These allowed the model to accommodate the carbon-electron temporal mismatch.
Eqn (1) of the previous section presents ordinary differential equations describing the system connectivity, while eqn (2) shows the functional representations. The model is more complicated than earlier models of glycolysis in L. lactis, as it addresses the pathway dynamics under anaerobic conditions. In contrast to aerobic conditions, the organism cannot easily recycle NAD^{+} under anaerobic conditions, and the balance between NAD^{+} and NADH therefore changes dynamically. These changes must be expected to affect the dynamics of glycolysis and are therefore considered important for the model.
Fig. 3B visualizes the dynamics of fluxes of accumulation (v_{1}) and consumption (v_{2}) of FBP relative to one another for the representative case of Experiment 3 with 80 mM of initial glucose. The horizontal top and bottom axes for the two curves are different and color-coded with blue (glucose) and black (FBP). In the first phase of the experiment, glucose is abundant and FBP accumulates (shown with a curved green upward arrow) while glucose is being consumed at a more or less constant rate (shown with the horizontal green arrow indicating that glucose is consumed at the same time as FBP accumulates). The intersection of the curves represents a temporarily steady state for FBP with a concentration of about 50 mM. At this state, accumulation and consumption of FBP are in balance: v_{1} = v_{2}, while other processes in the system are not. As the experiment proceeds, FBP remains constant, whereas glucose is being consumed from about 50 mM down to about 10 mM (color-coded with straight orange arrow). Since the value of v_{1} is essentially constant, FBP remains at its temporary plateau during this period. As soon as the glucose concentration drops below about 10 mM, FBP becomes depleted while glucose is being used up as the two red arrows indicate. Panel (C) shows the same three phases on a concentration vs. time plot of glucose and FBP for further clarification.
The concepts outlined above were implemented in a simple model. We start by representing the flux terms with Michaelis–Menten rate laws. In order to obtain the same peak for different amounts of glucose, v_{1} needs to be saturated and equal to V_{max1} for glucose with a value between about 10 to 80 mM for Experiment 1 and between about 5 to 40 mM for Experiment 2. These numerical settings permit the temporarily steady state and extended peak for FBP if v_{2} has a high V_{max2} and a low K_{m2}. In fact, these values must be such that for the peak FBP concentration (about 50 mM), v_{2} equals v_{1} for the appropriate glucose concentrations.
These quantitative considerations can be converted into constraints for the model parameters. Specifically, the equation v_{2}(FBP = 50) = v_{1}(5 ≤ glucose ≤ 50) mandates the following: First, K_{m1} must be such that v_{1} is saturated for high glucose values (above about 5 mM); this ensures a similar FBP peak for all pertinent glucose concentrations in the three experiments. Second, the FBP accumulation and consumption fluxes v_{1} and v_{2} need to intersect. This requires K_{m2} > K_{m1} and V_{max2} > V_{max1}, and the values must be set such that v_{1} = v_{2} when FBP is at the observed plateau of about 50 mM. A solution to these constraints is:
The relationships between fluxes are not dependent on the choice of the Michaelis–Menten framework. To reproduce similar conditions for power-law representations, one must simply require a low kinetic order and rate constant for glucose that results in the correct dynamics for v_{1}. The result is a similar flux value for different glucose concentrations between 5 and 50 mM. Furthermore, v_{2} needs to have a higher kinetic order so that the two flux curves intersect. Rate constants are calculated such that v_{1} = v_{2} for the peak FBP concentration of about 50 mM. Similar to the Michaelis–Menten case, the kinetic order and rate constant for v_{1} need to be set such that they satisfy the speed for glucose depletion. Results for a representative, feasible parameter set with power-law functions are illustrated with the fluxes in Fig. 4C. Fig. 4D shows the corresponding glucose and FBP concentrations for the three experiments with different initial glucose concentrations.
The question arises of whether the above considerations are affected by the existence of G6P as an intermediate between glucose and FBP. The short answer is no. Unfortunately, the [1-^{13}C] NMR experiments did not permit measurements of G6P for the three experiments. However, it is likely that the G6P dynamics is similar to the FBP dynamics, although with a lower peak value or a very low saturation threshold as it is observed for Experiment 1 with initial glucose of 20 mM.
Supposing that v_{3} is a function of FBP only, we find: (i) at times when FBP is constant, v_{3} is constant as well; and (ii) FBP being constant requires that v_{2} = v_{3}. (i) and (ii) are possible in the following two scenarios only when G6P and FBP are constant simultaneously, or when v_{2} is saturated for the amount of G6P during that time period, and therefore has a very low K_{m}. Furthermore, no matter whether G6P is constant or v_{2} has a low K_{m} (high affinity) for G6P, v_{1} = v_{2} = v_{3} must hold at the FBP peak, and v_{1} must thus be saturated for glucose values between about 5 and 50 mM as reasoned before. Thus, the earlier conclusions regarding the FBP dynamics hold, whether or not G6P is explicitly presented.
(3) |
(4) |
X0_{80} = [80; 0.1; 4; 14.8; 0.1; 0.1; 5.74; 0.1; 0.1; 8.815; 0.1; 0.1]; |
X0_{40} = [40; 0.1; 4; 14.8; 0.1; 0.1; 5.74; 0.1; 0.1; 8.815; 0.1; 0.1]; |
X0_{20} = [20; 0.1; 4; 14.8; 0.1; 0.1; 5.74; 0.1; 0.1; 8.815; 0.1; 0.1]; |
M_{B−20} = 13.92 mg protein per ml |
M_{B−40} = 17.11 mg protein per ml |
M_{B−80} = 19.53 mg protein per ml |
T_{20} = 0.7124 min |
T_{40} = 0.7450 min |
T_{80} = 1.6661 min |
These parameter settings lead to good fits of the data, which are displayed in panels (A–C) of Fig. 5. The capacity of a single parameter set to capture different conditions is important, because it significantly increases the predictive power of the model. Furthermore, since this parameter set was derived from DFE, its extrapolation reliability is increased, because the risk of error compensation among flux terms within the same or in different equations is greatly reduced.^{7}
While the simultaneous fitting with one parameter set is important, one should consider that, in reality, natural systems obviously exhibit a certain degree of variability. We therefore allowed the parameters in the common set, which fits the three experiments simultaneously, to vary slightly in order to account for this variability among the different cell populations. The results are shown in panels (D–F) of Fig. 5. These model instantiations used the following parameter sets for the 20, 40, and 80 mM experiment. A parameter-by-parameter comparison demonstrates how close the three sets are.
(5) |
(6) |
(7) |
For the first time, a single model fits all available metabolic time courses reasonably well with the same parameter set and also captures the counterintuitive FBP dynamics across the three experiments. Because this model reflects three independent datasets, one might expect that it has a higher extrapolation potential than earlier models that were based on single datasets. The companion paper^{1} discusses key aspects of the control of the pathway, and explains in detail how the organism terminates glycolysis at the end of glucose availability such that it is ready to take up now glucose as soon as substrate becomes available again.
This journal is © The Royal Society of Chemistry 2016 |