Monte-Carlo Modeling of the Central Carbon Metabolism of Lactococcus lactis: Insights into Metabolic Regulation

Metabolic pathways are complex dynamic systems whose response to perturbations and environmental challenges are governed by multiple interdependencies between enzyme properties, reactions rates, and substrate levels. Understanding the dynamics arising from such a network can be greatly enhanced by the construction of a computational model that embodies the properties of the respective system. Such models aim to incorporate mechanistic details of cellular interactions to mimic the temporal behavior of the biochemical reaction system and usually require substantial knowledge of kinetic parameters to allow meaningful conclusions. Several approaches have been suggested to overcome the severe data requirements of kinetic modeling, including the use of approximative kinetics and Monte-Carlo sampling of reaction parameters. In this work, we employ a probabilistic approach to study the response of a complex metabolic system, the central metabolism of the lactic acid bacterium Lactococcus lactis, subject to perturbations and brief periods of starvation. Supplementing existing methodologies, we show that it is possible to acquire a detailed understanding of the control properties of a corresponding metabolic pathway model that is directly based on experimental observations. In particular, we delineate the role of enzymatic regulation to maintain metabolic stability and metabolic recovery after periods of starvation. It is shown that the feedforward activation of the pyruvate kinase by fructose-1,6-bisphosphate qualitatively alters the bifurcation structure of the corresponding pathway model, indicating a crucial role of enzymatic regulation to prevent metabolic collapse for low external concentrations of glucose. We argue that similar probabilistic methodologies will help our understanding of dynamic properties of small-, medium- and large-scale metabolic networks models.


The Kinetic Rate Equations
In the following, we provide a comprehensive list of rate equations used within our analysis. We note that kinetic parameters with identical name across different kinetic equations (such as V max or K eq ) are meant to represent different parameters. In total, the system of differential equations encompasse 77 unknown Michaelis-Menten constants. We note that our choice of rate equations takes inhibition between substrates and products into account.  shows different behaviours of the systems depending on the specific set of sampled parameter values and the presence or lack of regulation. The ATP concentration during starvation can assume a wider variety of values in the presence of regulation and can also show oscillations. The average value of ATP concentration during starvation is 0.53mM for regulated systems versus 0.11mM for non-regulated systems. A common feature for both regulated and non-regulated system is the presence of a peak in ATP concentration at the time in which the starvation start (minute 1.0). The presence of this peak can be explained by the fact that, once external glucose has been deprived, the first part of glycolysis slows down for lack of its substrate resulting in a slower consumption of ATP by PFK. The downstream part of glycolysis might have a delay in responding to this stress, hence there will be a time-span during which ATP is produced faster than it is consumed.

Selected Time Courses of Metabolite Concentrations
As detailed in the main text, both in presence and lack of regulation we registered cases in which the system was characterised by bistability. For example, after restoring the pre-starvation level of external glucose, the 0.7% of the regulated systems and the 1.2% of the unregulated systems reached a new steady-state concentration of APT which was: (a) different from the original value within the 10% of tollerance and (b) sensibly different from zero (> 0.5mM). This phenomenon was in most cases mirrored by the hysteresys curves of such bistable systems. However, we also detected cases in which there was a discrepancy between the time course and what the hysteresys curve seemed to suggest. Fig. S3 shows two of such occurrences for the regulated systems. [ATP] (mM) Figure S2: Bistability of some of the regulated systems. The ATP concentration at steady-state after the recovering time is different from the pre-starvation value (within the 10% of tolerance) and from zero (> 0.5mM).
We explored this phenomenon further by generating hysteresis curves where we gradually changed the concentration of external glucose as external parameter.   S3 shows such hysteresis curves for the same models referred to in Fig. S2. In the case shown in Fig. S2-a and Fig. S3-a, we note that before starvation the concentration of ATP is on the lower side of the hysteresis curve, where the external glucose concentration of 20 mM corresponds to an ATP concentration of about 4 mM. After starvation the ATP concentration increases very steeply, reaching a peak of over 9 mM. From there the system relaxes to the new steady-state with the concentration of ATP decreasing until it reaches the upper part of the hysteresis curve, where the restored concentration of external glucose (20 mM) corresponds to a higher level of ATP (= 8.2 mM).
In addition to the results presented within the main text, we evaluated the control properties of a second metabolic state, with concentrations and fluxes given in Table S4. The second state was obtained using an explicit kintic model where all Michaelis-Menten parameters are set to values that are identical to the respective concentrations. Using this model, the value of external glucose was lowered to 2mM and the resulting steady state was recorded. In contrast to the metabolic state evalauted within the main text, the second state is therefor not directly based on experimentally acquired data. Nonetheless, analysis of the state allows for some conclusions with respect to typical control properties.

Reaction step Flux (mM/min) Species
Conc (  The probabilistic control analysis of the regulated and unregulated system was performed on the second metabolic state characterised by a low concentration of external glucose (see previous section for the complete list of concentrations and fluxes). For this second metabolic state we again observed that the presence of regulation was associated with an increased stability. In particular almost 98% of the regulated systems resulted to be stable versus the 84% of the non-regulted case.
We also notice that the control of the PTS on the glycolytic flux is positive in the second metabolic state for the regulated case . Fig. S5 shows the probabilistic sign distribution of flux control coefficients for regulated and non-regulated systems in the second metabolic state. Figure S5: Probabilistic sign distribution of the flux control coefficients for non-regulated (left panel) and regulated (right panel) systems in the second metabolic state. As in Fig.4 of the main text, the shade of each entry represents the percentage of the corresponding control coefficients that are positive. Figures S6A and S6B show the probabilistic distributions of the flux control coefficients for regulated and non-regulated systems in the second metabolic state. In this case the differences between presence and absence of regulation are more pronounces than in the metabolic state referred to in the main text.
In particular in the presence of regulation the distributions tend to be much slimmer, showing how the regulatory mechanisms substantially contribute in confining the value of the control properties of the system. This suggest that the action of the regulation is even more pronounced in a state with low glucose. Figure S6A: Probabilistic distribution of the flux control coefficients for the non-regulated system in the second metabolic state. As in Fig.3 of the main text, the distributions refer to the scaled flux control coefficients corresponding to the pathway map of L. lactis central metabolism. Columns denote changes in enzyme levels and rows the respective relative change in flux. Figure S6B: Probabilistic distribution of the flux control coefficients for the regulated system in the second metabolic state. As in Fig.3 of the main text, the distributions refer to the scaled flux control coefficients corresponding to the pathway map of L. lactis central metabolism. Columns denote changes in enzyme levels and rows the respective relative change in flux. Figure S7A: Convergence of selected control coefficients as a function of sample size. Shown is the median of 12 selected flux control coefficient C P i as a function of the number N of samples used for its estimation, with N = 1, . . . , 2 · 10 4 . We note that the scale on the x-axis is logarithmic. For sample sizes larger than a few hundred samples, no appreciable variation is observed. The figure corresponds to the control coeeficients shown in Figure 3 of the main text (full system, including regulation). Figure S7B: Convergence of control coefficients as a function of sample size. The figure corresponds to the flux control coeeficients shown in Figure 3 of the main text (full system, including regulation). Shown is the convergence of the median of all control coefficients as a function of sample size. The scaling of the axes is identical to those of Figure S7A. All axes span the interval [−1, 1] on the y-axis and N = 1, . . . , 2 · 10 4 on the x-axis. The latter is scaled logarithmically (See Figure S7A for labels). The flux control coefficient of changing the ATPase on five fluxes (PDH, ACALDH, ADH, PTA, ACK) is larger than one. These values are not shown, convergence properties are similar to those shown. Figure S8: Time-courses of intermediate metabolites. The median of the concentration of FBP, ATP, and PEP following a withdrawal of external glucose at t = 1min in models that do not posses regulatory interactions but are nonetheless able to recover after a period of starvation. The figure corresponds to Figure 11 of the main text. In contrast to Figure 11, all intermediate concentrations, including PEP, drop to low values after withdrawal of glucose. With the restoration of external glucose at t = 10min, PEP and FBP slowly rise while it takes until t ≈ 25 min before the concentration of ATP is restored to pre-starvation levels. The figure clearly shows that while fine-tuning of parameters might likewise, at least for a small subset of parameters, result in recovery after starvation, mechanism is very different from the recovery induced by feedforward activation.